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INTRODUCTION 


FVSOLVR  is  a two-dimensional  finite  difference  field  solving  algo- 
rithm for  the  calculation  of  electric  field  and  potential  distribution. 

The  nodal  solutions  from  FVSOLVR  can  be  used  to  plot  electric  field 
lines  and/or  equipotential  lines  using  the  CALCOMP  plotter  and  the 
companion  program  FVPLOT. 

These  programs  were  developed  as  part  of  a Navy  research  effort 
concerning  insulators  employed  in  large  antenna  arrays.  These  antenna 
insulators  are  subjected  to  severe  electrical  field  conditions  which 
occasionally  result  in  arcing,  corona,  and  heating  or  burning  of  the 
insulators.  In  order  to  determine  causes  of  insulator  failures  and 
evaluate  alternative  designs  it  is  helpful  to  know  the  electric  field 
distribution  associated  with  each  insulator. 

The  general  capabilities  of  the  programs,  the  theoretical  foundation, 
and  a discussion  of  possible  applications  can  be  found  in  the  Civil 
Engineering  Laboratory  Technical  Note  N-1502,  "Computer  Technique  for 
Calculation  of  Potential  in  Multidielectric  Media,"  by  Kwang-Ta  Huang, 
Brian  R.  Milner,  and  Andrew  W.  McClaine.  The  intention  of  this  manual 
is  to  provide  a user  of  FVSOLVR  and  FVPLOT  with  the  detailed  information 
required  in  order  to  implement  problem  solutions  with  these  programs. 

The  programs  described  herein  are  written  for  use  with  a Fortran  IV 
compiler.  The  program  versions  listed  in  Appendix  B have  been  used 
successfully  on  the  CDC  7600  computer  at  the  Lawrence  Berkeley  Laboratory, 
Berkeley,  California.  Users  are  cautioned  that  minor  modifications  may 
be  necessary  if  these  programs  are  used  with  other  systems. 


DESCRIPTION  OF  PROGRAM  FVSOLVR 

FVSOLVR  solves  problems  which  are  geometric  in  nature.  The  geometry 
of  the  solution  region  can  be  described  using  either  a cartesian  coordi- 
nate system  or  a cylindrical  coordinate  system.  The  solution  region  in 
the  cartesian  system  is  a rectangle  with  a side  parallel  to  the  x-axis. 

In  Figure  1,  the  x-axis  can  be  either  of  the  two  axes  shown  and  the 
solution  region  may  generally  include  parts  of  any  of  the  quadrants. 

(Some  versions  of  FVSOLVR  which  are  designed  to  take  advantage  of  the 
symmetry  of  certain  problems  require  that  the  rectangular  solution 
region  be  entirely  within  the  first  quadrant.)  The  rectangle  describes 
a plane  in  space,  perpendicular  to  the  z-axis,  with  z coordinate  z . 

The  condition  required  for  this  problem  to  be  two-dimensional  is  tRat, 
as  z varies,  the  field  solution  does  not  vary.  That  is,  any  similar 
rectangle  translated  in  the  plus  or  minus  z direction  will  have  the  same 
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solution.  See  Figure  2.  An  example  of  such  a problem  might  be  a capaci- 
tor plate  within  a housing,  as  shown  in  the  cross-sectional  view  of 
Figure  3.  The  z-axis  is  perpendicular  to  the  plane  of  the  paper.  The 
geometry  in  any  plane  parallel  to  the  paper  must  have  the  same  configu- 
ration. Idealizations  usually  must  be  made  in  cases  such  as  this  since 
the  capacitor  plate  and  housing  do  not  have  infinitely  long  z dimensions. 
But  if  d is  the  largest  gap  dimension  and  the  actual  z dimension  of  the 
capacitor  is  greater  than  20d,  then  experience  dictates  that  the  two- 
dimensional  solution  is  usually  a good  approximation,  until  one  gets 
closer  than  20d  to  a boundary  in  the  z direction. 

The  cylindrical  coordinate  system,  illustrated  in  Figure  4,  can 
also  be  used.  The  solution  region  is  either  a cylinder  of  radius  r = r j , 
or  a cylindrical  section  from  r = r^  to  r = rj.  In  either  case,  the 
solution  region  has  a height,  h,  which  is  measured  in  the  z direction. 

The  z-axis  must  be  an  axis  of  symmetry;  that  is,  the  geometry  and  solution 
are  two-dimensional  and  independent  of  0.  The  geometry  can  vary  with  r 
and  z and  the  function  to  be  solved  is  0 = 0(r,z),  with  30/96  = 0. 

The  program  incorporates  the  useful  characteristics  of  variable 
grid  spacing  and  variable  dielectric  properties,  thus  allowing  multi- 
dielectric medium  problems  to  be  solved  with  economy  of  computer  time. 

The  program  also  has  the  capability  of  solving  anisotropic  media  problems. 
Greater  detail  about  the  inner  workings  of  the  computer  program  and 
relevant  electromagnetic  theory  can  be  found  in  CEL  Technical  Note 
N-1502. 

Input  data  for  the  programs  are  developed  and  entered  in  the  same 
way  for  both  coordinate  systems.  When  using  the  cylindrical  system,  the 
z-axis  must  be  an  axis  of  symmetry.  Figure  5 illustrates  an  example  of 
the  solution  region  in  the  cylindrical  case.  The  known  circular  boundary 
in  Figure  5 will  be  explained  later.  The  r-axis  may  or  may  not  have  a 
plane  of  symmetry  passing  through  it  perpendicular  to  the  z-axis.  In 
the  cartesian  system,  the  y-z  plane  and/or  the  x-z  plane  may  or  may  not 
be  planes  of  symmetry.  Input  data  include  grid  coordinates,  dielectric 
values  (which  can  vary  throughout  the  grid  system),  dielectric  material 
boundaries,  and  known  potentials  or  electric  field  conditions.  Develop- 
ment of  program  input  data  will  be  shown  by  example,  using  the  cylindrical 
coordinate  system  and  solving  an  electric  potential  distribution  problem. 


PROBLEM  SOLUTION  USING  FVSOLVR 
Statement  of  Problem 

A potential  field  solution  for  the  region  surrounding  two  parallel 
corona  rings  is  required.  In  this  problem,  a potential  is  applied 
across  the  air  gap  separating  the  corona  rings,  which  are  metal  toroids. 
See  Figure  6.  Each  toroid  has  an  outside  diameter  of  20  inches  and  an 
inside  diameter  of  10  inches.  Each  toroid  is  suspended  by  a pipe  attached 
to  a thin  sheet  of  metal  welded  to  the  center  of  the  toroid.  In  this 
case,  the  z-axis  coincides  with  the  center  line  of  the  pipes  and  is 
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perpendicular  to  the  planes  of  the  toroids.  As  an  example  of  an  appli- 
cation, in  Appendix  A of  this  manual,  the  FVSOLVR  program  is  used  to 
solve  this  particular  problem.  The  top  ring  and  pipe  consists  of  con- 
ducting material  at  a potential  of  +100  volts,  and  the  bott''ii  ring  and 
pipe  linewise,  except  at  a potential  of  -100  volts.  Not  only  does  this 
problem  have  rotational  symmetry  about  the  z-axis,  but  it  also  has 
symmetry  across  a plane  perpendicular  to  the  z-axis.  The  plane  defined 
by  (z=0,r)  is  a mirror  symmetry  plane.  Thus,  solving  the  upper  right 
hand  quadrant  will  give  solutions  for  the  other  three  quadrants  as  well. 
This  solution  region  is  shown  in  Figure  7.  Significant  computational 
advantage  is  gained  when  solving  just  one  quadrant.  Only  one-fourth  as 
many  equations  need  be  solved;  therefore,  computer  costs  are  greatly 
reduced.  Problems  such  as  this,  which  have  mirror  symmetry  from  top  to 
bottom,  will  have  an  equipotential  line  about  which  this  mirror  reflection 
will  occur  in  the  two-dimensional  solution.  The  potential  of  this 
symmetry  line  is  always  half  the  sum  of  the  upper  and  lower  potentials. 
Since  the  upper  potential  is  +100  volts  and  the  lower  potential  is  -100 
volts,  the  voltage  at  the  center  (z=0)  line  is  0 volts. 

The  voltages  in  any  problem  can  always  be  normalized  so  that  this 
central  line  has  a voltage  value  of  zero  volts  simply  by  adding  or 
subtracting  a constant  value  from  all  known  voltages.  In  general,  the 
relative  voltage  levels  can  be  shifted  by  any  constant  value.  The 
FVSOLVR  program  requires  that  all  voltages  be  non-negative.  If  it  is 
desired  to  solve  a problem  involving  negative  voltages,  all  voltages 
must  be  shifted  to  eliminate  the  negative  values.  After  the  solution  is 
obtained,  the  values  can  be  shifted  back  to  the  original  levels  to 
obtain  the  true  solution.  The  shape  of  the  flux  and  equipotential  lines 
will  not  be  affected  by  this  kind  of  manipulation. 

Development  of  Data 

The  generation  of  the  input  data  is  performed  manually,  and  usually 
requires  the  preparation  of  a good  drawing  of  the  solution  region.  Note 
that  use  of  graph  paper  will  facilitate  the  preparation  of  the  drawing. 
Grid  lines  are  constructed  on  the  drawing,  keeping  in  mind  that  any  type 
of  boundary  must  consist  of  straight  line  segments  beginning  and  ending 
on  nodes  of  grid  lines.  The  circle  which  represents  the  corona  ring  in 
Figure  7 will  in  fact  be  a series  of  straight  line  segments  approximating 
the  circle. 

At  this  time,  we  will  examine  boundary  conditions,  which  are  known 
at  r=0  (viz.  3V/3r=0)  and  at  z=0  (viz.  V=0).  If  the  grid  region  does 
not  include  the  axis  of  symmetry  but  begins  at  some  radial  value  greater 
than  zero,  the  boundary  conditions  for  that  edge  of  the  grid  region  must 
be  given  explicitly.  (A  means  for  inputting  this  data  is  described  in 
the  section  which  discusses  the  input  variables  NPOTS,  POT,  NLINES, 
etc.)  But  as  r and  z go  to  infinity,  boundary  conditions  are  not  known. 

A means  for  treating  such  problems  will  be  presented  here.  The  theoreti- 
cal basis  for  this  approach  can  be  found  in  CEL  Technical  Note  N'1502. 
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The  boundary  conditions  at  large  values  of  r and  2 can  be  approxi- 
mated to  sufficient  accuracy  by  considering  that  at  large  distances  the 
corona  ring  gap  appears  like  a point  gap.  The  point  gap  arrangement  is 
a simple  geometry  which  has  a well-known  field  distribution.  At  large 
distances  the  electric  flux  lines  around  a point  gap  form  near-perfect 
circles.  Such  an  electric  flux  line  can  be  used  as  the  necessary  boundary 
condition  at  large  distances.  Thus,  at  some  distance,  which  experience 
dictates  to  be  10  to  20  times  the  distance  of  the  gap  between  the  rings, 
a boundary,  assumed  to  be  a circle,  is  introduced.  This  circular  boundary 
is  sh'^wn  in  Figure  8.  Inside  the  circle,  the  medium  is  air  with  the 
relative  dielectric  value  taken  as  unity.  In  the  region  outside  the 
circle,  the  dielectric  value  is  taken  to  be  zero.  Mathematically,  this 
means  that  the  region  outside  the  circular  boundary  does  not  affect  the 
region  within  the  circular  boundary.  At  least  two  grid  lines  are  needed 
outside  the  circular  boundary,  one  in  the  r direction  and  the  other  in 
the  z direction.  These  two  grid  lines  coincide  with  the  edges  of  the 
solution  region  shown  by  dashed  lines  in  Figure  8.  The  circular  boundary 
and  zero  dielectric  obviate  the  need  for  explicit  boundary  conditions  at 
the  edges  of  the  solution  region.  (Again,  it  might  be  noted  that  this 
circular  boundary  is  a series  of  straight  line  approximations.) 

Variable  grid  spacing  is  used  in  the  program.  This  is  extremely 
valuable  in  that  proper  use  of  such  a scheme  can  reduce  the  number  of 
equations  which  must  be  solved.  Experience  is  needed  to  predict  where 
fine  grid  spacing  is  needed  and  where  coarse  grid  spacing  will  suffice. 

A "rule  of  thumb"  is  that  regions  in  the  neighborhood  of  the  gap,  where 
small  radii  occur,  tend  to  have  the  highest  voltage  gradients  and  there- 
fore need  more  grid  points  for  good  approximations.  Regions  far  from 
the  gap  tend  to  need  fewer  grid  lines  for  good  approximations.  Experience 
with  the  program  and  actual  experimental  results  will  aid  in  arriving  at 
optimal  grid  designs. 

Boundary  lines  between  different  dielectric  regions  are  constructed 
within  the  grid  system.  These  are  composed  of  line  segments  constructed 
either  horizontally,  vertically,  or  diagonally  between  adjacent  grid 
points.  Regions  of  known  potential  can  be  assigiied  any  convenient 
dielectric  value.  This  freedom  is  allowed  because,  for  regions  of  known 
potential,  the  program  does  not  solve  equations  involving  dielectric 
values.  However,  for  completeness  of  the  data  arrays,  all  regions  must 
be  assigned  a dielectric  value.  It  is  most  convenient  to  assign  to  a 
region  of  known  potential  a dielectric  constant  equal  to  that  of  an 
adjacent  region.  In  such  a case,  there  does  not  exist  a dielectric 
boundary  between  the  two  specified  regions  since  the  dielectric  constants 
are  taken  to  be  the  same. 

If  a problem  includes  a conducting  medium  which  is  of  unknown 
potential,  the  region  occupied  by  such  a medium  must  be  assigned  an 
appropriate  dielectric  constant  value.  Satisfactory  results  can  usually 
be  obtained  by  assigning  to  the  conducting  region  a dielectric  constant 
which  is  very  much  larger  than  (and  nearly  infinite  with  respect  to)  the 
relative  dielectric  constants  of  the  surrounding  insulating  media. 
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Since  computers  cannot  represent  infinitely  large  numbers,  there  is  an 
upper  bound  on  the  value  which  may  be  used  to  represent  the  relative 
dielectric  constant  of  the  region.  If  the  value  chosen  is  too  large,  it 
could  result  in  over  or  underflow  problems  with  the  computer  during 
execution  of  the  program.  If  the  value  chosen  is  not  large  enough,  it 
will  not  appear  "infinite"  with  respect  to  the  other  relative  dielectric 
constants  in  the  solution  region  of  the  problem.  The  symptom  of  this 
malady  will  show  up  in  the  potential  distribution  solution,  where  the 
whole  of  the  conducting  medium  will  not  be  at  an  equipotential . The 
optimum  dielectric  constant  value  for  a conducting  medium  of  unknown 
potential  will  be  machine  dependent  and  the  user  should  determine  a 
suitable  value  by  trial  and  error  method.  Satisfactory  results  have 
been  achieved  with^a  CDC  7600  computer  using  a relative  dielectric  value 
of  the  order  of  10  . 

The  data  deck  format  will  now  be  explained.  Throughout  this  manual, 
in  the  discussions  pertaining  to  computer  program  variables  such  as 
HX(I),  HY(J),  NX,  NY,  NEPR,  NEPZ,  BDR(I,J),  BDZ  (I,J),  etc.,  any  reference 
to  X holds  true  for  r atjd  vice-versa.  The  same  relationship  exists 
between  y and  z.  The  information  for  the  data  deck  can  be  obtained 
almost  totally  from  the  drawing  of  the  problem. 

The  first  data  card  has  the  following  information: 

NX,  NY,  NEPR,  NEPZ,  IFSKIP,  IFBDY 
Format  614 

NX  is  the  number  of  axial  grid  lines. 

NY  is  the  number  of  radial  grid  lines. 

NEPR  is  the  maximum  number  of  dielectric  boundaries  that  might  be 
encountered  in  any  given  row.  A row  is  defined  here  as  a slot  between 
two  adjacent  radial  grid  lines.  See  Figure  9.  Note  that  there  must 
always  be  a dielectric  boundary  which  coincides  with  the  axial  grid  line 
of  maximum  radial  value.  For  the  grid  in  Figure  9,  the  proper  value  of 
NEPR  is  four  (4),  since  there  are  four  dielectric  boundaries  in  the  Ith 
row  and  no  other  row  has  more  than  four. 

NEPZ  is  the  maximum  number  of  dielectric  boundaries  that  might  be 
encountered  in  any  given  column.  A column  is  considered  to  be  a slot 
between  two  adjacent  axial  grid  lines.  See  Figure  10.  There  must 
always  be  a dielectric  boundary  which  coincides  with  the  radial  grid 
line  of  maximum  axial  value.  For  the  grid  in  Figure  10,  the  proper 
value  of  NEPZ  is  ten  (10).  This  value  is  determined  by  the  Jth  column, 
which  has  more  dielectric  boundaries  than  any  other  column. 

IFSKIP  is  an  indicator  used  to  test  input  data  (to  be  explained 
later).  Briefly,  if  IFSKIP^O,  the  program  merely  reads  the  input  data 
and  sets  up  the  boundary  conditions  so  that  the  user  may  verify  that  the 
problem  has  been  programmed  correctly.  If  IFSKIP=0,  the  program  executes 
the  solving  algorithm. 
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IFBDY  is  an  indicator  which  identifies  the  type  of  boundary  condi- 
tions which  are  to  be  used  for  the  top  (radial)  and  right-hand  (axial) 
grid  boundaries  shown  as  dashed  lines  in  Figure  8.  IFBDY=0  indicates 
that  the  boundary  conditions  are  to  be  treated  implicitly,  as  for  the 
circular  boundary  in  the  example  of  Figure  8.  This  means  that  it  is  not 
necessary  to  provide  the  potentials  along  the  top  and  right-hand  bounda- 
ries explicitly  as  input  data  (i.e.,  these  boundary  potentials  can  be 
ignored).  IFBDYj^O  indicates  that  potentials  along  the  top  and  right-hand 
boundaries  are  given  as  input  data.  This  is  usually  done  if  the  solution 
region  happens  to  be  a fine  mesh  in  a section  of  a previously  solved 
large,  rough  mesh  grid,  so  that  the  boundary  potentials  can  be  interpo- 
lated from  the  large  grid  solution.  (Further  details  are  given  later.) 

The  next  set  of  information  read  is  HX(1)  for  1=1  through  NX. 

HX(I)  are  the  r coordinates  of  the  axial  grid  lines.  The  HX(1)  must  be 
ordered  in  increasing  magnitude  and  must  all  be  greater  than  or  equal  to 
zero.  HX(1)  is  always  the  distance  from  the  axis  (line)  of  symmetry  to 
the  1th  axial  grid  line.  The  format  for  each  punched  card  having  HX(I) 
data  is  10F8.0.  When  a decimal  point  is  used,  a maximum  of  3 digits  to 
the  right  of  the  decimal  point  is  taken;  digits  beyond  the  third  decimal 
place  are  rounded  off.  Units  of  measure  should  be  scaled  to  conform 
with  this  format. 

The  next  set  of  information  read  is  HY(J)  for  J = 1 through  NY. 

The  HY(J)  are  the  z coordinates  of  the  radial  grid  lines.  HY(J)  must  be 
ordered  in  increasing  magnitude  and  can  be  greater  than,  less  than,  or 
equal  to  zero.  The  format  for  HY(J)  data  cards  is  also  10F8.0.  Again, 
only  the  first  three  digits  to  the  right  of  the  decimal  point  are  taken 
and  additional  digits  to  the  right  are  rounded  off. 

HX(I)  always  represents  distances  measured  in  the  radial  direction 
for  cylindrical  problems.  Axisymmetric  problems  must  always  be  oriented 
in  this  manner.  For  problems  represented  in  the  cartesian  coordinate 
system  where  no  grid  boundary  coincides  with  a symmetry  plane,  the 
problem  can  be  oriented  in  any  manner  with  respect  to  HX(1)  and  HY(J). 
However,  the  program  execution  will  be  faster  and  less  costly  if  NX  is 
less  than  NY. 

The  next  information  that  is  read  is  BDR(1,J)  for  J = 1 through 
NEPR,  for  each  I from  1=1  through  NY-1.  This  information  indicates 
where  dielectric  boundaries  exist  in  the  radial  direction.  The  program 
"searches"  the  rows  (radial  slots)  between  adjacent  radial  grid  lines 
for  dielectric  boundaries.  The  first  index  (1)  is  the  row  index,  which 
is  the  integer  value  of  a given  radial  slot,  starting  with  1 at  the 
bottom  of  the  region.  The  second  index  (J)  is  the  dielectric  index 
which  indicates  the  Jth  boundary  encountered  in  the  Ith  row  when  moving 
outward  from  the  axis  of  symmetry. 

The  assignment  of  the  values  of  BDR(I,J)  is  explained  with  respect 
to  the  example  shown  in  Figure  11,  where  NEPR=3  and  NY-1=4.  If  a boundary 
exists  axially  (parallel  to  the  axis  of  symmetry),  BDR(I,J)  is  given  a 
value  of  two  times  the  index  if  the  axial  grid  line  on  which  the  boundary 
occurs.  If  the  boundary  line  is  diagonal,  it  is  assigned  a value  equal 


to  the  sum  of  the  indices  of  the  two  axial  gr''d  lines  between  which  the 
diagonal  boundary  occurs.  This  will  always  be  an  odd  integer,  whereas 
all  axial  boundaries  will  always  have  even  BDR  values.  If  a given  row 
has  fewer  than  NEPR  boundaries,  dummy  boundaries  must  be  placed  at  the 
last  axial  grid  line  to  satisfy  the  computer  program.  Therefore,  in 
this  example,  NEPR=3  requires  that  a dummy  boundary  be  placed  in  row  3 
at  grid  line  6 to  complete  the  BDR(1,J)  array.  If,  for  example,  NEPR 
were  to  equal  four  and  NY-1  were  to  equal  fifty,  four  dielectric  bounda- 
ries would  have  to  exist  for  each  one  of  the  fifty  rows.  It  is  quite 
possible  that,  in  many  rows,  several  boundaries  will  be  "stacked"  at  the 
last  grid  line  in  order  to  satisfy  the  program. 

Another  example  of  the  assignment  of  the  BDR(I,J)  values  is  illus- 
trated in  Figure  12.  For  this  case,  NEPR=4.  The  format  required  for 
the  input  of  the  BDR(1,J)  values  is  2014.  (Thus,  in  the  present  example, 
all  of  the  BDR(1,J)  data  would  go  on  one  card.)  In  general,  NY-1  groups 
of  NEPR  integers  are  used  to  describe  the  boundaries  that  might  be 
encountered  in  the  radial  direction. 

The  next  information  read  is  a series  of  numbers,  KO,  KN,  KK,  and 
EPSR(K)  for  K = 1 through  KK.  These  numbers  are  used  to  assign  dielectric 
values  to  the  regions  identified  by  the  radial  boundary  search  just 
described.  The  information  is  governed  by  the  order  of  appearance  of 
the  dielectrics.  For  example,  in  Figure  12,  for  rows  one,  two,  three 
and  four,  respectively,  the  dielectric  values,  in  the  order  of  appearance 
are; 


Row 

1: 

1.0 

10.0 

1.0 

Row 

2; 

20.0 

1.0 

10.0 

1.0 

Row 

3; 

20.0 

1.0 

10.0 

1 .0 

Row 

4: 

20.0 

1.0 

A given  order  of  appearance  is  referred  to  as  the  dielectric  sequence 
for  the  particular  slot.  KO  is  the  index  of  the  first  row  in  which  a 
new  squence  of  dielectrics  is  encountered.  KN  is  the  index  of  the  last 
row  having  that  same  sequence,  with  the  provision  that  if  KN-KO  is 
greater  than  1,  all  rows  between  the  KOth  and  KNth  rows  also  have  the 
same  dielectric  sequence.  KK  is  the  number  of  dielectrics  in  the  sequence. 

The  first  card  to  be  read  for  this  series  of  numbers  contains  KO, 

KN,  KK  for  the  first  dielectric  sequence.  The  format  is  314.  This  is 
followed  (on  the  next  data  card)  by  the  values  EPSR(K),  for  K = 1 through 
KK,  of  the  dielectric  constants  in  their  order  of  appearance.  The 
format  for  this  data  is  16F5.0.  This  process  is  repeated  for  successive 
dielectric  sequences  appearing  in  the  solution  region  (until  KN=NY-1). 

For  the  example  of  Figure  12,  this  part  of  the  data  is  as  follows: 

K0,KN,KK  1 1 3 

EPSR(K)  1.0  10.0  1.0 
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K0,KN,KK 

2 

3 

4 

EPSR(K) 

20.0 

1.0 

10.0  1.0 

K0,KN,KK 

4 

4 

2 

EPSR(K) 

20.0 

1.0 

Since  the  sequence  of  dielectrics  in  this  example  is  identical  for 
Rows  2 and  3,  the  values  for  these  two  rows  can  be  read  in  together. 

This  "shorthand"  is  quite  useful  for  most  problems,  where  many  successive 
rows  will  contain  the  same  sequence  of  dielectrics. 

A similar  procedure  is  used  for  input  of  boundaries  in  the  axial 
direction  search.  The  BDZ(I,J)  for  J = 1 through  NEPZ,  for  each  I from 
1=1  through  NX-1  are  read.  The  rules  which  govern  the  BDZ  data  are 
analogous  to  the  rules  for  the  BDR  data.  The  program  "searches"  the 
columns  (axial  slots)  between  adjacent  axial  grid  lines  for  dielectric 
boundaries.  BDZ(1,J)  represents  the  Jth  boundary  in  the  1th  column. 

There  are  NX-1  groups  of  NEPZ  boundaries  read  in  format  2014. 

Another  set  of  KO's,  KN's,  and  KK's  are  then  put  into  the  data  set 
with  their  accompanying  EPSZ(K)  values.  The  EPSZ  is  the  dielectric 
constant  value  in  the  radial  direction.  This  value  need  not  be  the  same 
as  EPSR  at  the  same  point,  which  represents  the  dielectric  constant  in 
the  axial  direction.  This  allows  anisotropic  medium  problems  to  be 
solved. 

Boundaries  which  are  radial  do  not  appear  in  the  BDR(I,J)  array. 
Boundaries  which  are  axial  do  not  appear  in  the  BDZ(I,J)  array.  Only 
diagonal  boundaries  appear  in  both  arrays.  This  fact  is  useful  in 
checking  data  accuracy  using  the  boundary  check  program  (to  be  discussed 
in  a subsequent  section). 

The  quantities  NGX,  NGY,  VG,  and  GH,  which  are  explained  below,  are 
read  only  if  IFBDY^O.  The  IFBDY^O  option  is  for  solving  the  potential 
distributions  in  a finer  grid  region  (i.e.,  with  a finer  mesh)  in  the 
lower  left-hand  corner  of  a previously  solved  coarse  mesh. 

If  IFBDY^O  the  next  data  card  contains  NGX,  NGY  is  format  214.  NGX 
is  the  number  of  "coarse"  axial  grid  lines  in  the  "fine"  grid  solution 
region,  including  the  axial  boundary  of  the  fine  grid  region  (which  must 
also  be  a "coarse”  grid  line)  and  the  r=0  axis.  The  "coarse"  grid  lines 
are  those  that  were  used  in  the  initial  solution  of  the  complete  problem 
(i.e.,  with  IFBDY=0)  that  was  run  before  the  IFBDY^O  option  is  to  be 
executed.  NGY  is  defined  in  an  analogous  manner  for  the  coarse  radial 
grid  lines. 

Next,  the  values  GH(I)  for  1=1  through  NGX  are  read.  These 
values  are  the  coordinates  of  the  r intercepts  of  the  coarse  axial  grid 
lines.  The  format  for  GH(I)  is  8F10.3. 

The  voltages  VG(I),  for  I = 1 through  NGX  are  read  next.  These  are 
the  node  voltages  along  the  radial  coarse  grid  line  which  coincides  with 
the  top,  radial  grid  line  of  the  fine  grid  solution  region.  These 
voltages  are  obtained  from  the  previously  solved  coarse  grid  solution. 

The  format  for  VG(1)  is  8F10.3. 
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Ill  a like  wav,  the  coordinates  are  read  for  the  z intercepts  of  the 
coarse  radial  grid  lines,  viz.  GH(J),  for  J = 1 through  NGY , format 
8F10.3.  Likewise,  the  voltages  along  the  right  side  axial  boundary  of 
the  fine  grid,  VG(J),  for  J = I through  NGY  are  ready  in  format  8F10.3. 

NGX,  NGY,  GH(1),  VG(1),  GH(J),  VG(J)  are  the  only  optional  data  for 
the  program.  The  following  information  is  required  in  all  cases,  irre- 
spective of  the  value  of  IFBDY. 

NPOTS  is  read  next,  in  format  14.  NPOTS  is  the  number  of  known 
potentials  in  the  problem.  For  example,  if  part  of  the  region  has  a 
voltage  of  +100  volts,  and  elsewhere  there  is  a grounded  region,  with  no 
other  known  potentials  present,  then  NP0TS=2.  The  known  potentials  must 
always  be  greater  than  or  equal  to  zero;  no  negative  potentials  are 
allowed.  If  negative  potentials  exist,  the  absolute  value  of  the  minimum 
potential  must  be  added  to  all  known  potentials  in  the  problem  so  that 
there  are  no  negative  potentials  remaining. 

Once  NPOTS  is  read,  the  known  potentials  must  be  assigned  to  those 
nodes  at  which  the  potentials  are  known.  This  is  done  in  a loop  performed 
a number  of  times  equal  to  NPOTS.  The  first  information  needed  is  POT, 
format  F8.0,  which  is  the  potential  assigned  to  a region  (or  several 
regions)  of  known  potential.  The  "region"  may  also  be  a single  point 
(i.e.,  a single  node).  Next,  NLINES,  format  14,  is  read.  NLINES  is  the 
number  of  rectangles  needed  to  describe  all  regions  of  the  particular 
potential  POT.  These  rectangles  need  not  be  connected.  An  example  of 
how  to  calculate  NLINES  will  be  given  with  the  aid  of  Figure  13.  The 
fourteen  dark  dots  are  points  of  known  potential.  In  this  example,  it 
takes  two  rectangles  to  describe  the  region  of  known  potential.  Note 
that  one  of  the  rectangles  is  "collapsed"  to  a line  segment.  A point 
would  be  considered  as  a rectangle  that  has  "collapsed"  to  a point. 

Thus  NLINES=2  in  the  example  of  Figure  13. 

Next,  the  quantities  LI,  L2 , Ml,  and  M2  are  read  in  format  414.  LI 
and  L2  are  the  indices  of  the  axial  grid  lines  that  form  "edges"  of  a 
rectangle  of  potential  POT.  LI  must  be  less  than  or  equal  to  L2 . 

Likewise,  Ml  and  M2  are  the  indices  of  the  radial  grid  lines  that  form 
the  remaining  "edges"  of  the  rectangle.  Ml  must  be  less  than  or  equal 
to  M2.  So  for  the  example  of  Figure  13,  two  groups  of  four  numbers  each 
are  used  to  describe  the  location  of  the  two  rectangles.  These  numbers, 
read  into  the  computer  on  two  cards  following  the  value  of  NLINES,  would 
be  1 -2 , 1 -2 , J- i , J followed  by  1 - 1 , 1 +2 , J , J+2 . As  a further  example,  the 
values  of  1 1 , L2 , Ml,  and  M2  for  a single  point  of  known  potential  where 
the  22nd  axial  line  and  the  44th  radial  line  intersect  would  be: 

22  22  44  44.  Each  set  of  four  values  LI,  L2,  Ml,  and  M2  must  be  on  a 

separate  data  card. 

The  above  procedure  is  repeated  for  each  region  (which  may  or  may 
not  be  a "connected"  region)  of  distinct  potential  POT,  with  the  order 
of  input  data  being:  POT,  NLINES,  and  NLINES  groups  of  LI,  L2 , Ml,  M2. 

This  concludes  the  input  data  development  for  FVSOLVF;  a summary  of 
the  order  in  which  the  input  data  cards  are  read  in  the  FVSOLVR  program 
is  given  below. 


READ  NX , N Y , NE PR , NE PZ , I FSK I P , I FBDY 
Format  614 

READ  (HX(I).  1=1, NX) 

Format  10F8.0 

READ  (HY{J),  J=1,NY) 

Format  10F8.0 

READ  ((BDR(1,J),  J=l,  NEPR),  1=1, NY-1) 
Format  2014 

READ  K0,KN,KK 
Format  314 

READ  (EPSR(K),  K=1,KK) 

Format  16F5.0 

READ  ((BDZ(1,J),  J=1,NEPZ),  1=1, NX-1) 
Format  2014 

READ  K0,KN,KK 
F o rma  t 314 

READ  (EPSZ(K) , K=1 ,KK) 

Format  16F5.0 


The  following  five  READ  statements  are  executed  only 


READ  NGX,NGY 
Format  2l<* 

READ  (GH(l),  1=1 ,NGX) 
Format  8F10.3 

READ  (VG(1),  1=1, NGX) 
Format  8F10.3 

READ  (GH(J),  J=1,NGY) 
Formal  8F10.3 

READ  (VG(J),  J=1,NGY) 
Formal  8F10.3 


READ  NPOTS 
Format  14 

READ  POT 
Format  F8.0 


READ  NLINES 
Formal  14 

READ  LI, L2, Ml, M2 
Format  414 


] 


Repeated  NLINES  Times 


This  series  of  data 
cards  is  repeated, 
as  required  (see  the 
discussion  on  dielec- 
tric sequences). 


This  series  of  data 
cards  is  repeated, 
as  required  (see  the 
discussion  on  dielec- 
tric sequences). 

if  IFBDY^O. 


Repeated  NPOTS 
Times 
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Checking  Input  Data 


Once  the  data  deck  is  set  up,  FVSOLVR  can  be  executed,  but  it  is 
advisable  to  check  the  data  with  three  data  check  routines.  This  can 
save  time  and  money  by  catching  data  errors  before  the  solving  routine 
is  executed. 

The  first  routine  involves  a boundary  check  program  called  BDYCHK. 
This  program  plots  (using  Calcomp  plotting)  dielectric  boundaries  in 
such  a manner  that  a horizontal  boundary  appears  as  a horizontal  line 
segment,  at  the  given  location,  a vertical  boundary  line  appears  as  a 
vertical  line  segment,  and  a diagonal  boundary  appears  as  two  intersecting 
slashes  forming  an  X.  All  diagonal  boundaries  appear  in  both  the  BDR(I,J) 
data  and  the  BDZ(I,J)  data.  Therefore  the  appearance  of  only  one  slash 
of  the  X indicates  an  error  in  the  BDR(1,J)  and/or  the  BDZ(I,J)  data, 
since  one  slash  (\)  of  the  X comes  from  the  BDR(1,J)  data  and  the  other 
slash  (/)  of  the  X comes  from  the  BDZ(I,J)  data.  An  example  of  such  a 
plot  for  the  configuration  of  Figure  8 is  shown  in  Figure  14.  Numerical 
tables  are  also  printed  by  the  boundary  check  program  to  aid  location  of 
improper  data. 

The  data  deck  set  up  for  the  boundary  check  program  BDYCHK  is  as 
follows.  The  first  card  contains  NX,  NY,  NEPR,  and  NEPZ.  These  integer 
variables  are  defined  in  the  same  way  as  in  FVSOLVR.  The  second  card 
contains  10,  IN,  JO,  and  JN.  These  data  outline  the  region  to  be  examined 
by  the  boundary  check  program.  10  and  IN  are  the  indices  of  the  axial 
grid  lines  and  JO  and  JN  are  the  indices  of  the  radial  grid  lines  that 
form  the  borders  (or  edges)  of  the  region  which  the  boundary  check 
program  is  to  examine.  10  is  less  than  IN  and  JO  is  less  than  JN.  The 
HX(I)  cards  come  next,  followed  by  the  HY(J)  cards,  the  BDR(I,J)  cards 
and,  lastly,  the  BDZ(I,J)  cards.  The  cards  containing  the  HX(1),  HY(J), 
BDR(I,J),  and  the  BDZ(1,J)  data  are  the  exact  same  cards  as  those  prepared 
for  the  FVSOLVR  program.  In  summary,  the  deck  set  up  is  as  follows. 

READ  NX, NY, NEPR, NEPZ 
Format  414 

READ  10, IN, JO, JN 
Format  414 

READ  (HX(1),  1=1 ,NX) 

Format  10F8.0 

READ  (HY(J),  J=l,NY) 

Formal  10F8.0 

READ  ((BDR(I,J),  J=1,NEPR),  1=1, NY-1) 

Formal  2014 

READ  ((BDZ(I,J),  J=1,NEPZ),  1=1, NX-1) 

Form.ll  2014 


I 
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Once  the  boundary  data  are  verified,  a check  of  the  dielectric 
constant  information  is  performed.  The  given  dielectric  constants  must 
be  assigned  within  the  boundaries.  The  name  of  the  program  which  performs 
this  check  is  EPSCHK.  The  EPSCHK  program  uses  the  boundary  information, 
BDR(I,J),  BDZ(I,J),  and  the  dielectric  constant  assignment  information, 

KO,  KN,  KK,  and  EPSR(KK)  for  the  radial  direction  and  KO,  KN,  KK  and 
EPSZ(KK)  for  the  axial  direction,  in  order  to  generate  (as  printed 
output)  two  arrays  of  values  of  dielectric  constants.  One  array  contains 
the  values  of  dielectric  constants  encountered  in  the  radial  direction 
search  (rows),  and  the  other  array  contains  the  values  of  dielectric 
constants  encountered  in  the  axial  direction  search  (columns).  Two 
tables  are  printed  out,  one  for  each  array  of  dielectric  constants.  The 
two  arrays  are  compared  and  inconsistencies  between  them  are  noted  in 
the  program  output. 

The  data  deck  set  up  for  the  EPSCHK  program  is  as  follows.  The 
first  card  contains  NX,  NY,  NEPR,  and  NEPZ;  this  is  the  same  as  the  NX, 

NY,  NEPR,  NEPZ  card  used  in  the  boundary  check  program.  The  BDR(I,J) 
come  next,  followed  by  the  KO,  KN , KK , EPSR(KK)  cards,  set  up  in  the 
same  way  as  shown  in  the  FVSOLVR  program.  Finally,  the  BDZ(1,J)  cards 
and  accompanying  KO,  KN,  KK,  EPSZ(KK)  cards  are  read.  A summary  of  the 
order  of  reading  data  in  the  EPSCHK  program  is  as  follows. 


READ  NX, NY, NEPR, NEPZ 
Format  A 14 


READ  ((BDR(1,J),  J=1,NEPR),  1=1, NY-1) 
Format  2014 


READ  K0,KN,KK 
Format  314 

READ  (EPSR(K),  K=1,KK) 
Format  16F5.0 


Repeated,  as  required. 


READ  ((BDZ(1,J),  J=1,NEPZ),  1=1, NX-1) 
Format  2014 


READ  K0,KN,KK 
Format  314 

READ  (EPSZ(K),  K=1,KK) 
Format  16F5.0 


Repeated,  as  required. 


A third  data  check  is  used  to  verify  the  values  and  locations  of 
known  input  potentials.  This  check  involves  the  use  of  FVSOLVR  itself 
The  check  is  implemented  by  simply  using  the  data  deck  to  be  used  with 
the  actual  program,  except  with  IFSKIP^O  in  the  NX,  NY,  NEPR,  NEPZ, 
IFSKIP,  IFBDY  card.  This  results  in  FVSOLVR  skipping  execution  of  the 
solving  algorithm  and  printing  out  all  known  potentials  at  their  given 
grid  locations  and  "-1"  at  all  other  grid  locations 


i 
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Potentials  are  printed  by  a line  printer  as  well  as  written  on 
tape;  therefore  the  user  must  supply  the  tape.  Of  course,  if  the  user 
wishes  to  change  the  output  statements  in  the  program,  the  output  may  be 
obtained  in  any  form.  The  array  V(I,J)  contains  the  potential  solution 
at  the  nodes  having  coordinates  (HX(I ) ,HY(J) ) . 

Potential  listings  for  FVSOLVR  are  printed  out  in  the  following 
manner.  Blocks  of  50  values  of  potential  in  the  axial  direction  and  14 
values  in  the  radial  direction,  maximum,  are  printed  on  each  page  of 
output  from  the  computer.  The  table  is  started  from  the  lower  left  hand 
corner  of  the  solution  region  and  proceeds  to  the  right  producing  blocks 
of  14x50  potential  values.  This  is  illustrated  in  Figure  15.  If,  for 
example,  as  in  "page  3"  in  Figure  15,  fewer  than  14  values  of  potential 
are  needed  to  complete  the  printing  of  the  data  in  the  radial  direction, 
then  the  printout  will  fill  only  a portion  of  the  page.  A similar 
procedure  is  followed  in  the  axial  direction,  until  the  entire  solution 
region  is  printed. 
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USE  OF  THE  PLOT  PROGRAM  FVPLOT 

Use  of  the  plotter  program  FVPLOT  will  now  be  described.  The 
plotter  can  be  set  up  to  plot  equipotential  and/or  flux  lines.  However, 
as  a result  of  the  inherent  nature  of  the  plotting  program,  flux  lines 
can  be  plotted  only  for  problems  of  the  "point  gap”  type.  Examples  of 
"point  gap"  problems  are  shown  in  Figure  16.  Plotting  flux  lines  in 
configurations  other  than  the  "point  gap"  type  may  cause  difficulties. 
Equipotential  lines  can  always  be  plotted,  irrespective  of  whether  or 
not  the  problem  is  of  a "point  gap"  nature. 

Data  for  FVPLOT  may  be  input  from  cards  or  tape.  The  user  must  use 
the  appropriate  "read"  statements.  The  first  information  read  is  IVU. 

If  IVU=1,  the  program  plots  only  the  potential  lines.  If  IVU=2,  the 
program  plots  only  the  field  lines.  If  IVU=3,  the  program  plots  both 
potential  and  field  lines.  (If  flux  lines  are  to  be  plotted,  FVSOLVR 
must  be  executed  to  solve  for  the  necessary  field  values.  For  details, 
see  CEL  Technical  Note  N-1502.) 

Next,  a card  containing  10,  IN,  JO,  JN,  NQUAD,  KK,  VMIN,  and  SZ  is 
read.  10  and  IN  are  the  first  and  last  indices,  respectively,  of  the 
axial  grid  lines  that  form  edges  of  the  region  to  be  plotted.  JO  and  JN 
are  the  first  and  last  indices,  respectively,  of  the  radial  grid  lines 
that  form  edges  of  the  region  to  be  plotted.  NQUAD  indicates  whether 
only  the  solution  region  solved  for  by  FVSOLVR  is  to  be  plotted  or 
whether  "mirror  images"  of  this  solution  region  are  also  to  be  plotted, 
as  is  illustrated  in  Figure  17.  If  NQUAD=1 , only  the  solution  region 
obtained  from  FVSOLVR  is  plotted;  if  NQUAD=2 , one  mirror  reflection  is 
introduced;  if  NQUAD=4,  two  mirror  reflections  are  introduced.  This 
feature  enables  problems  of  varying  degrees  of  symmetry  to  be  solved  by 
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FVSOLVR  and  then  plotted  in  their  entirety  by  FVPLOT  while  not  necessarily 
inputting  the  full  geometrical  configuration  as  data  to  either  program. 

KK  is  the  number  of  equipotential  lines  to  be  plotted  "automatically"  by 
FVPLOT.  SZ  is  the  voltage  increment  between  successive  equipotential 
lines  and  VMIN  is  the  potential  value  which  is  incremented  by  SZ  to 
determine  the  first  equipotential  line  to  be  plotted.  Successive  incre- 
ments of  SZ  are  added  and  an  equipotenital  line  is  plotted  for  each  of 
the  KK  potential  values.  The  first  potential  value  plotted  is  VMIN  + SZ 
and  the  last  value  plotted  is  VMIN  + (KK  x SZ).  The  value  VMIN  is  not 
plotted  because  it  is  assumed  that  VMIN  usually  represents  an  electrode 
or  a boundary  of  the  solution  region.  Likewise,  any  known  voltages, 
such  as  electrodes,  should  not  be  plotted  "automatically"  as  just 
described.  These  voltages  should  be  given  as  separate  input  data  to 
FVGRPHR,  which  will  be  explained  later.  As  an  example,  if  two  electrodes 
of  known  potential  are  0.0  and  100.0  volts  and  equipotential  lines  are 
desired  at  10.0  volt  intervals,  then  KK=9,  VM1N=0.0,  and  SZ=10.0.  The 
0.0  and  100.0  volt  lines  are  not  plotted  "automatically"  and  are  not 
counted  when  determining  KK. 

NX  and  NY  are  read  next.  NX  is  the  number  of  axial  grid  lines  and 
NY  is  the  number  of  radial  grid  lines  used  in  the  FVSOLVR  solution. 
Following  this,  the  HX(1)  and  HY(J)  data  are  read  in  the  order  HX(1), 
HX(2),  HX(3),  ...,  HX(NX),  HY(1),  HY(2),  HYO),  ...  HY(NY).  Then  the 
voltages  V(I,J),  previously  calculated  (by  FVSOLVR)  at  the  grid  points, 
are  read  as  input  to  FVPLOT.  The  NX,  NY,  HX(I),  HY(J),  and  V(I,J)  data 
can  be  stored  on  tape  directly  from  the  FVSOLVR  solution  and  read  from 
the  same  tape  into  FVPLOT. 

If  IVU^l,  separate  data  must  be  input  for  the  field  (flux  lines) 
plot.  If  IVU=1,  the  separate  field  data  is  to  be  omitted.  The  field 
data  are  MX  and  MY  (which  are  exactly  analogous  to  the  number  of  grid 
lines,  NX  and  NY,  in  the  equipotential  case),  GX(I),  GY(J),  and  U(I,J) 
(which  are  exactly  analogous  to  HX(1),  HY(J),  and  V(1,J)  for  the  equipo- 
tential case).  This  information  can  also  be  stored  on  tape  from  the 
FVSOLVR  solution  and  read  from  the  same  tape  as  input  to  FVPLOT. 

The  last  information  inputted  is  that  of  various  reference  lines 
into  the  plot.  Such  lines  may  include  known  potential  lines  (as  previ- 
ously mentioned),  dielectric  boundaries,  or  any  other  lines  that  the 
user  wishes  to  plot.  Data  needed  to  plot  known  potential  or  reference 
lines  is  as  follows.  NLINES  is  the  number  of  lines  drawn  by  the  plotter 
to  generate  all  figures  and  lines.  In  determining  the  numerical  value 
of  NLINES,  a distinct  line  is  counted  each  time  the  pen  is  lifted  off 
the  paper.  NPTS,  which  is  read  next,  is  the  number  of  points  needed  to 
generate  a given  line.  For  each  line  that  is  to  be  plotted,  a data  card 
containing  a value  for  NPTS  is  followed  by  a card  (or  cards)  containing 
values  of  LX(N)  and  LY(N),  with  the  LX(N)  and  LY(N)  data  appearing  on 
the  cards  in  the  order  LX(1),  LY(1),  LX(2),  LY(2),  LX(3),  LY(3),  ..., 
LX(NPTS),  LY(NPTS).  NPTS  is  the  number  of  points  being  used  to  generate 
the  particular  line,  and  LX(N)  and  LY(N)  are  the  indices  of  the  axial 
and  radial  grid  lines,  respectively,  which  intersect  at  point  N.  The 
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plotter  draws  straight  lines  between  these  points.  If  a long  straight 
line  is  to  be  drawn,  only  the  end  points  of  the  line  need  be  given. 

Each  time  the  line  turns  a new  point  is  needed. 

The  following  is  a short  review  of  the  order  in  which  the  input 
data  cards  are  read  in  the  FVPLOT  program.  Each  READ  statement  by  which 
data  are  read  most  conveniently  from  tapes  is  indicated  by  an  asterisk. 

READ  IVU 
Format  14 

READ  IO,IN,JO,JN,NQUAD,KK,VMIN,SZ 
Format  614,  2F5.0 

*READ  NX, NY 
Format  214 

*READ  (HX(I),  1=1,  NX),  (HY(J),  J=1,NY) 

Format  8F10.3 

*R£AD  ((V(I,J),  1=1, NX),  J=1,NY) 

Format  8F10.3 

The  previous  five  READ  statements  are  executed  irrespective  of  the 
value  assigned  to  the  integer  variable  IVU.  The  following  three  READ 
statements  are  executed  only  if  IVU=2  or  IVU=3. 

*READ  MX, MY 
Format  214 

*READ  (GXd),  1=1, MX),  (GY(J),  J=1,MY) 

Format  8F10.3 

*READ  ((U(I,J),  1=1, MX),  J=1,MY) 

Format  8F10.3 

READ  NLINES 
Format  14 

READ  NPTS 
Format  14 

READ  (LX(N),  LY(N),  N=1,NPTS)  Repeat  NLINES  times 
Format  2014 
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QUADRANT  II 


QUADRANT  I 


Figure  1.  Labeling  of  quadrants  in  the  cartesian 
coordinate  system. 


QUADRANT  III 


QUADRANT  IV 


100  V 


Capacitor  plate  within  a housing.  The  geometry  of  the 
capacitor  and  the  field  solution  is  the  same  in  any  plane 
perpendicular  to  the  z— axis  (the  z— axis  is  perpendicular 
to  the  plane  of  the  paper) . 


Figure  3.  Example  of  a problem  described  In 
a cartesian  coordinate  system. 
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Dielectric  boundaries  are  shown  as  heavy  tines.  Hie  numbers  in  the 
j’.rtd  are  relative  dielei'.ric  constants.  That  is,  in  the  Ith  row, 
the  relative  dielectric  constant  chanRes  from  6.0  to  20.0  to  1.0  to 
20.0  in  passlnR  outward  from  the  axis  of  symmetry  in  the  radial 
direction.  Tliere  must  bt'  a dielectric  boundary  on  the  griii  line  of 
maximum  radial  value. 


fiRure  9.  Kxamples  of  dielectric  boundaries  encountered  in  the 
radial  direction. 
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Dielectric  boundaries  are  shown  as  heavy  lines.  Tlie  numbers  in 
the  grid  are  relative  dielectric  constants.  Tliat  is.  In  the  Jth 
column,  the  relative  dielectric  constant  changes  from  20.0  to  1.0 
to  20.0  to  ...1.0  to  20.0  to  1,0  In  going  from  bottom  to  top  in 
the  axial  direction.  There  must  be  a dielectric  boundary  on  the 
grid  line  of  maximum  axial  value. 


Figure  10.  Examples  of  dielectric  boundaries  encountered  in 
the  axial  direction. 
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axial  grid  line  indicies 
3 4 5 


The  heavy  lines  are  dielectric  boundaries;  and  the 
numbers  in  the  grid  are  dielectric  constant  values 


BDR(1,1)=4 
BDR(2, 1)=5 
BDR( 3, 1)=6 
BnR(4, 1)=5 


BDR(l ,2)=8 
BDR(2,2)=8 
BDR(3,2)=12 
BDR(4,2)=8 


BDR(l , 3)=12 
BDR(2,3)=12 
BDR(3,  3)  =12 
BDR(4,3)=12 


Figure  11.  Numbering  scheme  used  in  the  assignment  of  the  BDR 
(I,J)  values  where  radially  encountered  dielectric 
boundaries  are  specified. 
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'llie  heav'v  lines  are  dieleetrie  boundaries;  and  the 
numluTs  in  t ie  urid  are  dieleetrie  eonstant  values 


Bl)R(l,l)=6 
BDR( 2.1)=) 
BDR(), 1)=A 
BDR(4, 1)=) 


BI)R(  1 ,2)  = I0 
Bl)R(2,2)=b 
B1)R(  ).2)=6 
BDR(4.2)=14 


BDRd,  3)=14 
BnR(2. ))=10 
B1)R(3,))=10 
B1-)R(4,3)  = 14 


BDR(1 ,4)=14 
BDR(2,4)=14 
BDR(3,4)=14 
BDR(4,4)=14 


Figure  12.  Kxample  of  the  use  of  several  dummy  boundaries  at 
the  last  axial  grid  line  to  eomplete  the  BDRCI.J) 
array . 


"rectangles 


Figure  13.  A region  of  known  potential  "POT"  which  Is 
described  bv  two  "rectangles". 
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j In  this  example  NX=34  and  NY=118.  There  exists  correspondence 

between  the  position  of  data  in  the  printout  and  the  location 
‘ in  the  solution  region.  The  solution  region  is  indicated  by 

a dashed  line  in  this  figure.  In  this  particular  case  nine 
pages  of  computer  printed  output  are  required. 

• ‘ ! 

^ Figure  15.  Example  of  the  format  of  the  printout  of  field 

. values  from  FVSOLVR. 
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Figure  i6.  Kxamples  of  problems  of  a point  gap  nature 
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Figure  17.  Illustration  of  the  effect  of  the  NQUAD 
assignment  statement. 


NQUAD=1 

1 Quadrant  Plotted 


NQUAD=2 

2 Quadrants  Plotted 


NQUAD=4 

4 Quadrants  Plotted 


Appendix  A 
EXAMPLE  PROBLEM 
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An  example  of  a potential  distribution  problem  is  given  to  illustrate 
the  organization  of  the  input  data  for  FVSOLVR  and  FVPLOT.  This  simple 
problem,  which  is  described  in  the  text  and  illustrated  in  Figures  6,  7, 
and  8,  involves  only  two  dielectrics  and  two  given  potentials.  One  of  j 

the  dielectric  values  is  zero  and  is  used  only  to  satisfy  boundary  i 

conditions  near  the  top  and  right-hand  edges  of  the  grid.  Many  problems 
which  are  commonly  encountered  are  more  complicated  than  this  example, 
and  will  require  a more  elaborate  input  data  set.  However,  the  data 
format  will  be  Just  as  in  this  example. 

The  corona  ring  and  pipe  are  assigned  a potential  of  100  volts. 

The  radial  line  at  the  bottom  of  the  solution  region  (plane  of  symmetry) 
is  assigned  a potential  of  0 (zero)  volts.  The  region  inside  the  circular 
boundary  is  assigned  a relative  dielectric  constant  of  1.  This  includes 
the  regions  of  the  corona  ring  and  pipe  which  can  have  any  arbitrary 
dielectric  constant  since  the  potential  in  those  regions  is  known.  The 
value  of  1 is  chosen  for  convenience.  The  region  outside  the  circular 
boundary  is  assigned  a relative  dielectric  constant  of  0 (zero).  Thus 
it  is  not  necessary  to  assign  boundary  potentials  to  the  top  and  right- 
hand  edges  of  the  solution  region. 

The  grid  lines  in  the  solution  region  are  shown  in  Figures  A-1  and 
A-2.  Table  A-1  lists  the  input  data  for  FVSOLVR.  Table  A-2  lists  the 
input  data  for  FVPLOT  for  the  equipotential  line  plot  near  the  corona 
ring  as  shown  in  Figure  A-3.  Figure  A-A  is  an  equipotential  line  plot 
of  the  entire  solution  region. 
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The  grid  lines  near  the  corona  ring  are  too  densely  packed  to 
allow  good  resolution  in  this  figure.  Figure  A-2  Illustrates 
the  grid  lines  near  the  corona  ring. 


Grid  lines  for  the  entire  solution  region  for 
the  corona  ring  problem  shown  in  Figure  8. 
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TABLE  A-1.  INPUT  DATA  FOR  FVSOLVR 


NX,  NY,  nEP2,  IFSXIP,  IFtfor 
60  61  2 1 0 0 


HX{I) 


0.00 

0.50 

1 .00 

1.27 

1.55 

1 .6U 

2.51 

3.18 

3.80 

9.50 

6. as 

5.1/0 

5.051 

5.179 

5,«29 

5,776 

6.  lUU 

6.520 

6.735 

7.003 

7.500 

7.997 

6.265 

8.UH0 

8,656 

9.22U 

9.571 

9,621 

9,9U9 

10.00 

10. O'* 

10.26 

10. «9 

10.76 

11,09 

1 1 .60 

12.02 

15.23 

1U.08 

15.00 

17.20 

20.70 

27.50 

5«.60 

57,00 

71,10 

90,20 

113. OC 

155.70 

160,90 

189.50 

218.10 

256.30 

2UU.50 

252.50 

256.50 

260.50 

26U.56 

269.65 

265,*" 

HV(J1 


0.00 

9,00 

8.00 

15.00 

15.00 

19.00 

20,60 

21.12 

21,87 

22.55 

22.69 

23.21 

25.50 

25.75 

25.90 

29,0  0 

29.05 

29.12 

29.20 

29,90 

29,69 

25.10 

25.57 

26.110 

26.50 

27,00 

27.93 

27,90 

28.31 

28.60 

28.60 

26.66 

28.95 

29.00 

29.12 

29,38 

29.79 

30.28 

30.69 

31  .59 

32.50 

33.20 

39.00 

56,00 

98.63 

58.53 

60.93 

102.20 

119.95 

150.52 

165.29 

210.99 

227,62 

239.70 

299.18 

255.26 

258.60 

262.17 

263.60 

269.50 

265.00 


»UR( J,J} 


lie 

120 

118 

120 

116 

120 

118 

120 

117 

120 

116 

120 

116 

120 

116 

120 

116 

120 

116 

120 

116 

120 

116 

120 

116 

120 

116 

120 

116 

120 

116 

120 

116 

120 

116 

120 

116 

120 

116 

120 

116 

120 

116 

120 

116 

120 

116 

120 

1 16 

120 

116 

120 

116 

120 

1 16 

120 

1 16 

120 

116 

120 

116 

120 

116 

120 

116 

L’O 

116 

120 

116 

120 

1 1 6 

120 

116 

120 

lib 

120 

116 

120 

116 

120 

116 

120 

116 

120 

US 

120 

119 

1 20 

113 

120 

111 

120 

109 

120 

107 

120 

105 

120 

103 

120 

101 

120 

99 

120 

97 

120 

95 

120 

93 

120 

91 

120 

69 

120 

87 

120 

85 

120 

102 

120 

KO,  AN,  KK 


1 5«  2 


E8PR(kk) 


1.0  0,0 


38 


TABLE  A-1. 


(cont . ) 


K0»  KN,  HK 
kO  60  1 


E8P*>(MN) 

0.0 


eoz(i.j) 


120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

12o 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

120 

122 

110 

122 

117 

122 

1 15 

122 

115 

122 

111 

122 

109 

122 

U'7 

122 

105 

122 

103 

122 

101 

122 

99 

122 

97 

122 

95 

122 

95 

122 

91 

122 

67 

122 

1 1 

122 

122 

122 

KO(  »*««  MK 

i S6  2 


£SPZ(KK) 


1.0  0.0 


AO|  "N,  AH 
bP  4«  I 


ISP/IAA) 


0.0 


NPOT» 

i 
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TABLE  A-1.  (cont.) 

POT 

0.0 

nli^-es 

j 

1 

LJ.  Ll, 

1 60  1 1 


POT 


H 


100. 


NLI»<t8 

21 

Cl#  C2.  "I.  '*2 


1 

2 

2a 

61 

s 

S 

2S 

61 

H 

12 

2S 

2S 

11 

13 

29 

26 

1<< 

la 

23 

27 

IS 

IS 

22 

26 

16 

16 

21 

29 

17 

17 

20 

30 

IS 

la 

19 

31 

19 

19 

16 

32 

20 

20 

17 

33 

21 

21 

16 

3a 

22 

22 

17 

33 

21 

23 

16 

32 

2a 

2tt 

19 

31 

2S 

2S 

20 

30 

2b 

26 

21 

29 

27 

27 

22 

26 

2B 

26 

23 

27 

29 

29 

2a 

26 

30 

30 

2S 

2S 

TABLE  A-2.  INPUT  DATA  FOR  FVPLOT  FOR  FIGURE  A-3 


IVU 

I 


lOi  IN,  JO,  JN,  MJU*0,  KK,  V«IN,  SI 

i uu  1 as  t a 0.0  10. 0 


NX,  NV 

(OAT*  UN  TAPt  •'vSUUVN) 

HX(1),  HYIJ) 

(DATA  UN  TAPE  FPUN  FvSijLvR) 

V(1,J) 

(DAT*  UN  TAPE  FPJN  FvSdLVN) 

NLlNES 

J 

NPTS 

i 

CX(N),  l.Y(N) 

1 1 aa  1 
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TABLE  A-2.  (Cent.) 


NPT8 

5 

LX(N),  LY(N) 

1 2a  2 2U  2 25  3 2S  3 US 

NPTS 

30 


LX(N)»  LV(N) 


3 

25 

12 

25 

13 

20 

14 

27 

15 

21 

34 

22 

33 

23 

32 

2« 

31 

25 

2P 

2tt 

20 

23 

2T 

22 

26 

21 

25 

10 
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Appendix  B 


PROGRAM  LISTINGS 


The  Fortran  program  listings  for  FVSOLVR  and  FVPLOT  which  follow 
are  current  as  of  the  date  of  publication  of  this  report.  These  listings 
supersede  those  given  in  CEL  Technical  Note  N-1502.  Revisions  are  made 
periodically  as  deemed  necessary;  however  these  changes  are  usually 
concerned  with  minor  details.  Several  versions  of  FVSOLVR  and  FVPLOT 
exist.  For  sake  of  brevity  only  one  version  of  each  is  listed  here. 

The  version  of  FVSOLVR  given  here  is  applicable  to  geometries  which 
have  axial  symmetry.  Other  versions  apply  to  two-dimensional  cartesian 
coordinate  systems  which  may  or  may  not  have  planes  of  symmetry. 

FVSOLVR  uses  three  tapes  (or  disks)  during  execution.  Tape  33  and 
tape  3A  are  used  only  as  "scratchpads"  during  execution  of  the  algorithm. 
Upon  completion  of  the  program  execution  tape  10  contains  NX,  NY,  HX(NX) , 
HY(NY),  and  V(NX,NY),  the  last  array  being  the  nodal  solutions. 

The  version  of  FVPLOT  given  here  yields  Calcomp  plots  for  11-in. -wide 
paper.  The  plots  are  drawn  such  that  the  smallest  dimension  of  a rectan- 
gular plot  is  scaled  to  a width  of  10  inches.  Field  solution  data  from 
FVSOLVR  are  read  from  tape  10  (potential  solution  data)  and  tape  12 
(flux  solution  data). 

For  completeness,  the  optional  data-check  program  listings  are 
provided.  BDYCHK  aids  in  verification  of  placement  of  dielectric  bounda- 
ries. EPSCHK  aids  in  verification  of  assignment  of  dielectric  constants 
to  the  appropriate  regions. 
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PiT*  *p/>r,noP*n,p/,sv«/pPn*i’,i'/,B/?oo*ft,n/,v/2nnoo*"l.o/ 

POPmxT  ST*TF-r* TS 

1 POPPiTf toxp.ni 
j rn»“*Tf?Piu) 

X pnpfc*  . (i*.xs,n^ 

(I  Fnow,Tf«ri  P.X1 

S FnP-»T  f 1 M J , UM»  V ■ , Tu , « X , «h».y  *,iui 

t,  rroM*T(iM  ,*i  .XHMxf  ,Tx,xMi  *,*'S,x,a(*tt,6X,?MMx  f , rs,  b,fs,^)J 

7 FOPBXTnu  ,i1  , X«Hvf  , T J,  JM)  »,*-S.X,«(»u,6X,XMHvC,T?,3B)  a.FS.Sl) 

a FfiPMiTdMj  ,pit,  )M*,i  aPa.x^ 

p Pf|PM»T  f 1 nx,  1 B*  , 1 <iF*,  X ) 

IP  rnp»'»TfiPx,i»**,i«(»M.. 

1?  pnPMXT  f 1 xn , 1 butmti  4L  viiiiF  ■«sa.3) 

JX  PPPX*T  n PTS  «,T71 

IP  rnp-*T n X ,TMX.ii«mK«, TS1 
IS  F0O***Tf  1 x,SM»ipif5«,  TS» 
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100  Mcni.»MrnL»i 
lNKBNCnL*SO 

130  XlBKJ+t 

B^SKt  «TMK 

PBTNT  , fULwB. «»«*(•<  1 ,'<»'«  I .''?«sn) 
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300  FPS2(I.K)BFPS7(Kn,K) 

360  TF  (KK.EO.MPP7)  CnTO  BOO 
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C IMTEPPOLPTF  PTME  goto  POUMOPPT  PnTFMTHLS 
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00  BSO  Tal.MK 
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TF  (nsym.gt.O)  Goyn  5000 

TEPfl  40  4nO  5 4OO4YS 

on  0*00  j*i,Nv? 

5( JJaO.O 
on  0*0  0 1*1, Nn' 

0*00  4PTl,J1*0,0 

C41.L  54N0F0  N4TPTX  sni  VFP 
N*  1 

TF  fM«tt1.iF,fl.01  “a? 
r4LL  P4NSn(  (»0,NnPN,Nf.?1 
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OHTM  10 

OPTNT  9 , f MV f K ) , KrMl , N?1 

TA  OPTn  5^00 

TP  (M?,LT,A'yl  RfiTn  5100 

CnMTTMlif 

OTOP 

PMO 

OUPPnilTiME  0 AnjSni.  ( A , Ml  , N? ) 
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1.  SMiPT  91  nry  OF  fqijatipims 

MBsNBtl 
on  J?5  MSI, MW 
MMsMMOM 
9(M)s9f MM) 

5(MMls0,0 


( 
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o o c> 


JHIS  FAOE  IS  SE3I  ' 


<^n  I ?■;  Us)  ,u>i 

9,  oFtp  tjF«T  pmr<  nr  fqii4tip's  T^T^ 

tF  ( MIIMRI.  K.f  0 ■)  tFn,?0n,IFO 

OF*n  f^u)  (B  CO  . 4 f t , ►J)  . 4 f ?,  " 1 » 4 f “►*,  , WS^L  , MH) 

TF  (>.'0  9(1^,  ICICI, ?pn 

F,  WFniiCP  BLCriC  HF  FOIUTTI'mS 

pD  jnp  V*  1 , ^^' 

TF  C4CI,4')1  99S,^nn,??«; 

R(M)sBCNV4f',M) 

no  ^75 

TP  CAfl  ,'')1 
r»4(l  ,MW4C),|.T 
T s S ♦ I,  • 1 
TbP 

PO  ?F0  KKt,'“' 

I*  J*  1 

4(,I,T)s4CJ,T5-r»4(K,M 

BCT)sBCIT-4(L ,►  )»BCnT 

4(l.,»ij»c 

PO^TTnijf 

rnMTTMiF 

u,  nQTTF  PLpru  OF  PFnijrFp  Fi.ii4Ttn»'s  n>u  tapf  ? 

7F  (mimblk-mpi  %7«:,tt(lP,^7S 
C4I  L Bn'.BT^XCB,4Cl,t»,»(^,u“«i; 

nn  in  loP 

B4rF.S"PBTTTIiTTn‘i 

on  asn  p»i,4'4 

kJa^^4  t 

on 

I s M 4 K . I 

R(M)*S(^).4CF,^)*BCL'' 

M u * N * M KJ 
B(S,U7,BfW) 

TF  C*IB,FO.MU''BI  F,  AMn.v.RT.MCKPT  007(1  USP 

i.xsLy«t 

TF  fl*.CF.L««)  corn  aoS 
I,  VbLV«1 
I.  FbMC-I 
VCl  *,L7T»bCO 

COmTT'jiIF 

NBsNBal 

TF  (NB5  U7F,«;pn,U7S 

rOMTTMi.'E 

PAUL  »n«BT^^CB. 4C 1 , 1 1 ? ) 

on  Tn  OOP 

rnk'TTi'i'F 

BFTuBV 


FMO 
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THIS 


C 

Hl.mHOOT  IK'E 
OI^E^STPV 
oi^EMSTON  f 

PIAT4  KlUi/ft/ 

TF(^W,NF.01  Gr  T'^  G 
C*IL  nPE^^•s(T•^,«Tny,uoft,n) 

KPaM 
K*»N*“ 

S TPflMny.MF.t  1 r-n  to  io 

C >.OTTf 

1 

r*LL  KPIfMSf 

NwaM**  1 

CALL  »pIT-«(!;-\,a,ka,m«) 

GP  TP  ?n 

C bFaP 

to  CALL  OFAPMSf^,,*,<» 

CALL  PPAriMSt^-A.R.KB.k'Wl 
^wa^M« ] 

c nn^r 

?0  CP‘"tTMIIF 

9E  TUBW 
Fvr 

c 

SUPBnuTI^^E  »RT10fA,Mt 
m^EMSTOM  AfM) 

'■(OTTFtIO)  A 
OET  UPN; 

PnH 

C 

SU«BnuTI»J£  PBtV(A,B,S) 
ni-EMSTPW  P(V) 

PPTNT  11, A, P 

>1  FOP“ATf  iy,FP.'^,?w  ♦,1uFP.51 
OETUP'J 

FKin 

C 

FUBRnuTIMF  MPp I ? ( L V , J , FPS 1 ,WA»9B,ANEPS,A0FA,LF#R) 

C CALCi'LATFS  PTFI  FTTBIC-ARFa  PPUnuCT 

C PF  PPPITP'/TAL  P^ACF  nF  VPLIJPt  feLE“FNT 

c 

PT“EN8TPN  FPSdRQ.IOt 
TAiTEGFB  POonQP.IO) 
fP“PP^/MP9*BR/BnB,FP8 
t TF  (POPd  r. *t  1,?,1 

1 t»J»1 

FP81»EP8fLV,>»1 
TF  (Pr»P(LV,Jt-?*L*t  1,?,8 
? ,I»  I*1 

FP8?»FPSfL'',.») 

« A^FP8Bn.8*( (PpS1«EPS?)*P**?tFPS?*Rn**8>EPS1 ABAAa^) 

FP81 bFPS? 
rFTuP^ 

8 ANFP8BFP8t*APFA 

I •FTllPM 

* F«(P 

c 

8U»'»PUTT»/E  vFpT(L*.E*SI  A2EP8,LV,  V) 
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THIS  PA3E  IS  3£3I  aUALITI  PIUA>riA%B|^ 
^JU  OOFY  nOWi6HBD  K>  I»C  ^ 


! 


C rALCDLATFS  Oln  FCT9ir-»RE4  ^'Rl  niTT 
C nf  VFRT1C*L  F»rF  flF  volume  fcLtMENT 

c 

ni-EMSlON  FPSfOR,ift5 
TmTERER  Bn7(Rg,(n) 
roMMnK'/vFR4RR/pr)7,EPS 
K«1 

TF  (Bn7(i  1,?,-? 

1 K*«+1 

FPFl»EPSfL*.xi 
TF  (Rn7(t X,K )-?«L Vi 
P 

PPF?»FPSaF,X) 

47FPS»R*rFPS1*fv-71»FP5?*f7PM-Yn 
FPF1 *EPS? 

Rf  Tl'RN< 

^ 47FPS*FPF1 *R*f7PM-71 

return 

FNin 

C 

SuPonuTikE  FPrPEF  (M 

C r4l  ri'L'TFS  COFFirtfjTS  of  CpIC  SPI.Inf  uSEO  Tk  Intfrooi  4 TiOk 

C 

01 ‘'E MS. ON  f?ft,ii,Fv  fpnni,S(P>^n),«Miif?nni,T4U(?oo) 
rnMMnN/sP4P(;/v' . f*  , s 

N“ 1 UN* 1 

‘|KP«»J.p 

RMOfPIsn.fl 
T4ilf  Pisft.o 

on  1 

T“i »T-t 
TP1 »T»i 

Hl''JSVNfT).K»fT''11 

Mls*ki(TPl  l-xk  fT) 

TEMPsfwj«i/HTW(PH0tT1*?.)*,>, 

0M0( T*1 1 «.1 . /TP"0 

o«6,*ffF*)(TPn*F<.(Tn/MT-fF'-fn.F*i(t-p>/HT“n/MT 
1 T4ll(T*n«(r'-MT“t*T4liri1/Hp/Tt“P 

s(  n*p . 

Ff» 1»P. 

on  ^ Tst . ' M? 

TPxN-r 

? S(Tm«0«n(T«*i  1«4(  TM*I  1*T4Ii(  T'-*!  1 
RF  TUOV 
FNO 

c 

FUNCTrn-  $a|  7k, F fkj.xl 

r.  FV4LII4TFP  4P1  TNF  HINTTlnN 
C 

oi-Fn^tom  r?nft7,F.fpfi  ,7  >,«(?.  oi 
roMHPi\/sP4Br, /y».,fn,« 

TF  (v.r,F,«Nfp7  r.nin  1 

M 1 H M ( P 1 . F N f M 

SPl  tMF»FN(p*fy»y«ifp)*((F*  (?)-F«  f1  ))/mi««»«S(P1/6.) 

Rf  Tiidn 

1 TF  (F.LF.xkfNii  r,oTn  p 
>|M  1 SM«  1 

MU'.  I »«M  f 4 1 . «*.  f - t 1 

SPlIVE«FMf»,  )*rV«VNfNl)*rfFk(Nj-FNfN“1)WUN''i*k.4.41«S(>J>*P/*(,) 
RFTL'RN 
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OOiJ 


rOMinuB  I I'F  PLOTTFB 

’^I-EMSTnn'  v(ifli',?('n>,"*fi0«),'‘V(?nft),iPi/P(A0),G)(flft0),Gvf200), 

rn“MnN/0B»*'4Br,/<  Cbn('1.''t^^nn),»k■*x,yM4K,^^u4D^FTHSTX,nELT»X,FI98Tv 
/r>El.  T4V 
I FnPM4Tf BF1 
? rnPMjT  f?niu) 

^ F0PM4Tf6l«,?Fi;.n) 
u Fn9M4Tfia) 

s FnB'‘4T  »,  TU,h»,UM4iV  »,lu) 

«,  cnP«4T  ( IMfi  , T M,  T1  U,  T?U,  T1  u) 

7 FOBMaT  C7M01I  FBnM,Ft''.^,ttM  T'  , F 1 0 . 3 . fl  ><  , 6HV  F POM  , F ) Q . 3 . « « TO, 

BEAD  PnTF^TI4(,  r,BIPI  n4T4 
OEAP  (I  ,IVII 

BEAD  3 , TP,  T»  , JP,  J'J.  ^^l■4^,•  K . W-T'j,  S7 
'|PPT*S7 
OF 40  (tnT 
TALL  BPinfw»,>7T 
FAIL  «P1 P rPv , « V ) 

Pri  BO  J*  1 , f-v 
PB  CALL  OPi  n ( w f 1 , .1)  , N*  1 
TSCAKJ*) 

TF  (TVM.FQ.lT  r.'TTfi  ISO 
BEAP  Ftl'V  r.BTP  P4T4 
O£40  n?) 

FALL  ePl 0 f P* , T 
FAIL  BDIOfny.uvl 
pn  BO  .T*i  , “y 
00  FALL  BPI  0 f Ilf  1 , T)  , «*  T 
MM»| 
k|M«2 
KPAr  1 

Smtft  pata  Fro  off.toki  tp  of  blpttep 

IF  (TO.fP.n  PHTP  1?0 
pn  10  0 T s I r , T K, 

F.xriti-inisr.Vf  T) 
pn  1 oo  .t«i  . “V 
too  Ilf  T + 1-Tn,  Jlrlif  T,J) 
i?o  >4*«l>u*i.Tn 

TF  (JP.FO.n  F.nTP  ISO 
PO  lUO  Ts.jn,.!*' 

RTf  J*l-jn)»r,vf.T) 
pfi  lao  rn,*** 

1 UO  nf  T , .1*  1 -.TPlslif  T , T1 
ISO  TF  (Tn.FO.n  F.PTn  170 
PU  TfcO  I*in,Tt 
yi*f  T*l-Tr)»K»f  n 
pn  too  .T«i,*v 
lf.0  V(T71-Tr),  JTiVf  T,  IT 
170  MX3lM*1-Tn 
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o o o o r»  r> 


XBIS  PiSE  IS  BEST  ^UAL^  PIUCTICABU 
PBOB  COri  fUiij'ir'jtiED  XW 


^r  ( JP.F'^.i  1 RPTH  >nfi 
nn  loo  ,j»jn,,iv 

Dfi  i»>n  in  , Kit 
1 v(  T , j* t .jin av- r T,  n 

?oo 

s ,‘'y,»Y 

BBTN'T  , TP,  T>.,  jn,  Jki 

POTM  T ,HX  f 1 I .MVI  1 ) ,'-V(kiVl 

yYAXaNKfi^xI-w^d  1 

yM***MV(NVl»»-yf  1 1 
MX-l  ■•jy.1 
My**!  «Ny-i 

T»<TTI*LI7F  pi  PT  *».p  pc*.<  iYH 

C*l  L Pi  nts  (1i,o,PP) 
r»iL  Pint  fn,,.ni,,-si 
r*LL  PLnT(0.,r.S,-^) 

FIRSTX*m*( 1 I 
r j PSTy*MV ( 1 ) 

IF  (*«4x,GT.y«»*I  GOtn  ppn 

TF  (MQii* n ,F f,’,?.  n.  f ?. .x"* y ) .UT.yw**)  Gntn  ??e 
PEL  TiKay^FX/l ft. 

PFLTFyiOFLtAy 
Rntn  2«ft 

3?n  PFI T*y«y“»y/i  ft, 

PEI,T*X*OFUT*y 
TF  (kil5ll*n.Fn.?'  GPTO  PPft 
PUft  TF  (M(3U*P,RT,n  C.nTn  ?*ft 
TF  (y-iy.GT.y-iyl  Gntp  ?^n 
y*y«fy»*4y/y*<»»wio, 

CALL  AyiSfP., t''..?tH7.AyTS  AXJS  PF  SV«“E tpy , ♦^7, y 4 y , ft , , F TpSTy 
JP£I  TAX, 01 

TALL  AyTSfO.,ift.,*-Hp-»yTS,-ft,lil.,270,,FIQSTy,DFLT4y,ftl 
Gotn  ?oft 

?».ft  y4y«(y“Ay/y“4yi*i  0, 

TAIL  4yiS(ft.,ft.,P7H7»»yiS  AXTS  flF  ST-XF TPy , *27 , I ft , , PO , , F TP* Ty 

inEl.TAY.ft) 

TALL  »yrs(ft.,ft., ftwP.AX  TS, .ft. yAy,ft.,FrPSTy,ftFlTAy,ft) 

GHtp  PPO 

PPft  PFi  TAy»2. apfi  tax 
OFl TAysP,*PFI  TAX 
ppft  TF  (TV'I.FO.PI  n-OTP  ^ftS 

PB4W  EQuTPPTFftTlAl.  (.TftFS 

nn  Sftfl  K*i,yy 
7«v“ift*y*S7 

70ft  PAIL  PPAW  r 7, » y, vy , MX-I , Nyxl , py ,My , V, isr AM) 

TF  fTVIl.FO.n  70S, soft 

c 

C PALCULATF  AMP  PLPT  flux  I 1 ftt  L T 8 TP T S' IT T PM 

c 

70S  xyMi«My.i 
Myx) ixy-i 
T«1 

OAPVfl  , 1 1 
71ft  TSTAI 

BB*vn , n 
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THIS  PACtE  is  best  QUAXjITI 

jp,OM  oon  iTjeuJteHKD  ro  coo  — ^ 


Tf  (IIPnT.GF,P».4Mn.(lPOT,LF,MRJ  GOTO  ^?n 
P«aPn 
GOTP  J10 
52n  I^-al-l 

OmKY  f !“)♦  (UPPT-Pl ) *f  H¥ 

Many 
TSC*»Ja? 
raft  , 

P«aHX (KQi  } 

•5J0  KBRaKO»*l 

TF  (KOP.GT.x*)  GPTn  Rfio 
QBaHy(KBB) 

FRa(VCFOP.FHl.Vtl^Ba,*<H))/H 
t2B( (FP-F*5/f ) /?. 

*|aF4-P**?.*4? 

TF  (P.LT.OP)  rT'0C3*?*(B**?-«'***?)  + »l»(R.B*  > 

CaC*ri‘C 

TF  fP.LE.oo)  r,nrn  juft 

KBtawBR 

BAxBR 

r.nrn  jjn 

RUft  7*"(KP*,MH)*fP«B4)*(Ufl'BR,«>'1-H(«B»,“Hl)/CPR»B4) 

riLL  ,*■'¥“1  ,u*,Gy,ii,  TRr»NT 

JSft  OXaB 

rsuf-aft. 

^<>ft  *2a(  f FP«F»  Vf  UP-B41 ) /?. 

41 aF  4-B4*2,*4? 

CtFiCa4?«fBP»*?-By**?1*4l*(or'.yy) 

tf  (r-rsnM.cT«r)  jRft.jPn.jT" 

^7ft  rs'|*‘ac^u“♦rI^^ 

KB43KBP 
XBRaKB  4 * 1 

TF  (¥SB,GT.“yi  cnr.i  sftft 

P4aBP 

PlaP4 

OBaH* (KPO ) 

r4a(vfaP4,‘M>,vf«B4,‘'M))/H 
cntn  jfeo 

XBft  4ftaCSu*'»TC*4?*PX**?*4|*oy) 

r»I.L  Pjni  V f B,  B*  , OR,  4?,  4 i , 4ft  ) 

GHTO  UOfi 

Jpft  OaPR 

aOO  7B(i(aP4,»'H)  + ro-B4)»fiir«pa,“>')-u(KO4,»*rt))/f0a«B4) 
C4LL  DB4«  f z. >*,».¥, *-¥“1  ,My«l  , I**, r.y.ij,  TSC4N1 
TF  (Tsr4N,Ff;.«)  r,oTii  son 
TF  fP,K;F,OP)  ftO'Xn  jsn 
aP4aKPa 
BBRaKPa* 1 

TF  (KPa.GT.yyi  GnTO  Sftft 

PRaHv ( M BR ) 

P4aP 

GOTO  JRft 

nB4«  CLFrTPrrp  4»'ri  oTFLFftTPlC  FnnFTGii04  TTPiu 

Sftft  04  40  J I*JFS 

ftp  TOft  Lsl,*.LT*'FS  ,, 


JJHSP1.5E  IS  BEST  QUALITT  PRACTICABU 
jgiQH  UOi''i  lUSHkoHj©  XO 


. I 
r . 


i 


* 


9E*0  ^ ,NPTS 

PE40  ? , (I  Xf^)»L VfS?,N»1 ,nPTS) 

IF  (XN*x,GT,v»<*)0  r.rtTn  *,00 

TF  (NOUFO.fO.p.A^fn.  (?,*X“»X)  .1»T,  VH*X>  ROTH  600 

on  SttO  Kmtt^PTS 

l.xN^(^lsLxr^J)♦1»IO 

UYMMNPLvrNiti.jn 

v(M»hv(lvmk'X/i 

rjnTo  (51  n,sio,R?n,s?n)  ,NOii<iO 
510  X(N)»XM4x-MX((^*SMigl 
ROTH  5«0 

5?o  y(w)pVM*x.Hv(i  vMvtgl 
530  X(M)»X«*X*H* (t yvMMl 
5ao  rnMTiNiiF 

X(MPTS*1 5«FIP>TX 
XfK'PTS^JlPnFLTXX 
V(MPTStn*FTBpTV 
v(sPTS*?)«nFLT*r 
CALL  LINE  (Y.X.WPTS.I .0,0) 

RtlTO  C700,F«0,«150,550) 

550  on  560  ^*1,\pt.5 
560  X(M*?,*XM4*-x(m) 

rxLL  LINF  (V.v.MPTS.l  ,0,0) 

on  5T0  ^sl,^'PTS 

570  V(N)»?.*VM*X-V(N) 

C*I.L  LINE  (v,x,.jPTS,1  ,0,0) 

500  on  500  M«1,K/P7S 
500  xCKi)»?.*y*'4y-xfN) 

r*i  L LTNF  (y,y,»'PT8, 1 ,0,0) 

R07n  700 

600  on  6U0  ^IPI.NPTO 

LXMNN«LX (N) ♦ 1 »in 
i,yNNw»Ly(MUi.jo 
y{N)«My  (Ly^^N^) 

Rnm  (61  0,630, 6?0,6?0) 

610  X (M)  «HX  (L  XMI^^  ) 

Rorn  6ao 

6?o  y(M»y“4x-Hrf|  yMNKi) 

630  y ( •><) «x*‘»y4Hx f I y 'wvM ) 

6«0  rOMTTMUF 

y(NPTS*i IxFtouTx 
y(*PTS*?)»OFi  T*y 
y (I^PTS*!  )»FIB«TV 
y (6PT«*?)«r'Fi  T4y 
0*1. L LI»iF  (x,y,^PTS,i  ,0,0) 

ROTO  (70O,6«n^(,SO,65n)  ,fjRM*U 
650  on  660  1 , k Pt5 

660  y(Ms?,*yMxx>x(N) 

o*lL  (x,y,^0TS,  1 , 0,0) 

on  670  N»1,mOT5 

670  y(‘>)a?.*yM4x»y(\) 

PAIL  UT^F  (X,v,K'OT5,  1 ,0,0) 

600  on  650  NB 1 , ^PTF 

690  X('i)B?.»XM»i(.*fM) 

r*LL  UtMF  (x,y,NPT5,| ,0,0) 

700  rn».TTNHF 

r*i L pinT(ia.n,o.p,uo) 
oynp 

F^iO 


c 
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su'i'sniiTi^F.  tjniPfA.M 
ni^FMSTOk) 

PF40  (in)  * 

OFTLiPk. 

FVjPi 

c 

SliPPOijTIMF  f'P/ii-'  r7  . K , >L  » T T • Jj  , «V,  yv,  V,  ISC»^> 

C R«I''  ikin  PI  nTS  l,T».FS  I F Frjilii  V»LUF 

C 

''I“F.MSTn*''  yvfi''0'),vv/f?nnl,vn''0,?of') 

r(i>'‘inM/nP»i.*Pr,/»f*.no),yf*,n<M,»“4y,YyAy,'inu*0,FTPST**i'FLT*x,Fii*sTv, 
/f>Fi  T*v 
MTsn 
TysO 
T»7 
T4»IL 

nn  5 j«i  , J.I 

^AcJ 

43V ( T 4 , J 4 1 

''sVf  14,  J4*1  1 

TF  (4,fQ.ni  r.n  in  \ 

TF  C4.1F  .O..np.r'.LF,n.)  RdTl'  i 
TF  (4.F0.T1  1.? 

t T34«  f r)>4  ) * 1 , F.4i 

? call  ri4  (n,  4 , t,  y V , vv,  * , y , T4,  J4  • Ty  , mt  , ,jj,  VI 
TF  (MT.N'F.n)  r,r  Tn  io 
^ rn^-TTwiiF 
T43JL 

nn  6 l3t . IT 
tasil-i 

R8V(T4*1 , J4) 

43V(I4,.I4I 
TF  (R.FQ.AI  r,n  Tn  R 
TF  (4.LF.n..no.B.LF.0,1  RPTi.  h 
TF  (B.F0.T1  U,«. 

It  TeB*  f 4-P  ) * 1 , F-fr 

S CALL  AP ( 4 ,P, T , y V , y V, y , y , T4 , J4 . Ty , ^T, T T , VT 
TF  (NT.NF.ni  r.n  tr 
<1  rokiTTMilE 

TAlO 

on  10 

t43jL-J 
civd  ,.Ti*n 

RtVCl , J4) 

TF  (C.FO.RI  r.n  TO  10 
TF  (jA.FO.n  r.nTo  ? 

TF  (C.LF.0..nB.B.LF.0.1  ROTL'  10 
7 TF  (f  .LF.0,.rp.P.LT.O.l  Rf)TI’  10 
TF  (C.FQ.TI  F,0 
• T3r4.(R»Cl  *1  .F.f, 

0 CiiL  RC(P,r,T,yv,yv,x,y,T4,J4.iy,M,  Tj,vi 
TF  (MT.MF.o)  r.n  Tn  ?o 
10  rnMTTMUF 

f4«0 

on  la  r«i • rr 

TAlI 

riBVCTA.n 

r»v{T4*i , 1 ) 

TF  (O.FO.C)  c,n  Tn  1U 

TF  f TA.EQ.l  .»4O,tsC4M,F0.2)  G^Tn  il 

59 


9 


TF  (c,Le.o..no.n.LF.o,)  cnin  m 
11  TF  (C,LE.O..nP.O.tT.O.)  GtlTO  lu 
TF  (O.FO.TT  1?,13 
1?  T»OtfC-05«l ,F-h 

n CALL  COfr,n,T,*v,yv,x,v,T*,J*.ix,NT,i!,V) 
TF  (NT.MF.ni  r,n  to  sn 
la  roiyTiMiiE 

TF  (TSC*N.^f  .1  T cnm  ?? 

Tin  21  i«i,tl 

nn  2t  J»i , jj 

T4«I 

.UbJ 

C»V(T*1»J*I  ) 

F»V(T+1 , J) 

TF  (C.FO.fl)  r.f>  TD  Ifl 
TF  (.i.FQ.n  F.nrn  is 
TF  (C,LF.O..np.S.LF,0.1  ROTCI  lA 
is  TF  (C.LF.0..PP.H.LT.0.1  r.nTI)  lA 
TF  (C.eo.TT  U.17 


16  T*C*fS»C)»l ,F.6 


Fa 


mcntuaia 

I»r  i y JHJ*J^,riJSLl  TO  DOC 


TF  (T.FO.Vf  I*,  141)  '/(TA,J4)sVf  T4,JA1*1  .00f»l 

TF  (T.Fo.vfTA  + i , J4l ) vri44i »,)4)av(T4  + i , j4)*i ,nnni 
TF  ( T.FQ.  V M **1  » '44  1 1 ) V n 4*1  , JA+ 1 1 rV  r T4  + 1 , J4^  1 ) *1  ,(jnOl 
TF  tT.F-3.V(TA,J4*l))  VfT*,.I**nsW(TA,,I4  + 1)*1.n001 
43V( T4, J4) 

HsV  f T4  + 1 , ,J4) 

C * '/  ( T 4 ♦ 1 , ,1  A ♦ 1) 

) 

TF  (TA.FC.I  .4Kn.TSC4V.FfJ.?1  r.l  TO  ■?? 

TF  t.IA.F.O,  1 ) r.OTn  -44 

TF  ( A .1  F .n.  .f'B.H.LF.n.  .no,r.i  t .0.  .ns.o.i  F.n.  1 GOToao 
X?  TF  (4.1  T.o..ro.q.LF.n,.rD,r.|  t .o..n^.r).LT.o.i  GfiTn  so 
<;oTn  Xtt 

5X  TF  (4.LT.o..PQ.<4.i.T.n,.no.r.i,T,o..nB.n,LF.n,)  r,nTr.)  So 
5«  TF  (mt.fo.i)  r.r  to  ss 

TF  fk'T.FP.?)  r.n  Tfi 

TF  (MT.F^.X)  r-P  TO  ST 

TF  (MT.FP.U)  r,p  TO  5P 

RO  Tn  AO 

5S  PAU  n4(ri,A,T,VV,yV,K,V,T4,JA,TF,>jT,.IJ,V1 
TF  (Tf.RT.f'TF)  C.ptp  SQ 

PAIU  Bf('>,r,T,yy,VV,y,V,T4,JA,TX,M,.JJ,V) 

TF  (Ty.GT.MTV)  PfPTn  S<9 

P4I.L  4fifA,B,T,y\/,VV',y,V,TA,J4,I*,wT,II,V) 

GO  Tn  SQ 

X*1  pall  ABf  4 ,b,  T,y  V,  yv,  X,  V,  t a,  Ja,  tx  ,^T,  1 1 , V) 

TF  (TX.GT.HIVl  GHTO  SB 

TALL  Cn(C,o,T,*v,yY,y,y,T4,J4.Tx,MT,T!,V) 

TF  (IX.GT.KiTX)  GPTP  SB 

pall  Rr(P,r,T,xy,yv,x,y,TA,J4.TX,NT,Jj,y) 

GO  TO  SB 

JT  PALL  “P(B,P, T, yy.yv, X,y, T4, J4 . TX ,MT, JJ, V) 

TF  CTX.P.T.MIX)  Rpyp  SB 

PAI.L  n4(n,A,T,yy,yv,y,y,T4,J4.Tx,^T,JJ,v) 

TF  (Ty.GT.FJT*)  GHTn  SB 

pail  p0fr,0,T,yv,yv,y,y,i4,J4,Tx,MT,it,v) 

GO  Tn  so 

50  call  CPtP.o.T.yv.yv.y.y.TA.JA.Tx.MT.Tt.v) 

TF  (Tx.gt.mtx)  Rnsn  sb 

PALL  4«{A,B,T,yv,yv,y,y,TA,J4.Ty,^T,Il,v) 

TF  (IX.GT.^IX)  p.nTn  SB 

CALL  OA(P,4,T,yv,yv,y,y,TA,JAtTx,NT,.ti,v) 

50  TF  (TA.EO.  TST.AnP.  Ja.FO.  TSTJ  !jO,SI 
an  TF  ( y ( Ty  1 ,Fo,y  f Ty.n  1 gpto  So 
s»fvfly)«yfiy.l))/(xfiy)-y(iy»l)) 

SBtyf ry)»S*xf Tx) 

U1  0L«(yVf,IA*1)-yV(.lAl)/(XVfIA)»XV(TATn) 
0««(yVfJA*n«VV(JA))/(yvfIATlJ»XV(TA)) 
nOTHL»AOS(SL*1 ./S) 
n«TMB»4BS(S0T1 ./S) 
tF  (noTMO.LT.roTWL)  GDTn  a? 

SfiO«yV(J4).Sl  *XV(IAt1  ) 

*o«(san-sB)/fs-sL) 

RMTO  as 

(If  SftOavvf JA)»SB*yvrTA) 

yo«(5Pn-SB)/(s-SO) 

as  TXalVAl 

VdXTsVO 

v{ Tx)«s*y f txJtSH 

so  CALL  xPLnTV(y,y,7, Tx) 
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L, 


I 

i 


TF  (VMix.GT.VMdd  RrtTO  fcO 

IF  fNf3IIA0,FQ.?.*Hn,f?.*XM4K)^UT,V»<4X)  RnTD  hO 
nn  54)  i«i , Ty 

GOTO  fSt  ,55,S?,s?i  .nGiiapi 
51  *m=yM*y-yni 
GOTO  So 

5?  V( T )aVM*y-V f T 1 

55  »(T)*yMAy*yfl) 

50  rnMTNi.ife 

»(  I*  + nsFIPSTy 
y(Ty*2)snFI  T*y 
y(TX*l IsFTOSTv 
v(Ty+?IsnFi Ttv 
CALL  Llf'Fiy.y.ix.i  ,n,nl 
Gnin  fA,G,s«.sc,ss)  .v^itAo 

55  -yn  56  T»1 , Ty 

56  y f T ) s?.*y**A*-»  r n 

CAI  L LT*iF(v,y,Ty,i  ,n,p1 
no  57  1*1,1* 

57  y(T)*?,*y'<Ay-vrii 

C*IL  LT^iFCy.y,!*,  1 ,o,n1 

GO  50  T»l  , T * 

50  y(i)s?.*yKAy-y(n 

CAI  L Lr»-F(y.*,  I*.  i.G.G) 

GOTO  60 

60  GOTO  (65,61  ,61 ,61  ) ,'.n|,An 

61  on  6?  1*1,1* 
y(T)**“A**y(Ti 

6?  TF  (Nf.'liAn,FG.O  ) Vf  T)*?.«VMA».Y(n 
65  y(TX*n*FT0STy 
y(T*A?)«nFLTAy 
VCIy4.n«F10STy 
v(Ty*?)*nELTAv 
CALL  LI''lF(y,V,TX,1  .0,01 

GOTO  (69,67,60,601  ,K)GllAn 
60  nn  65  T*i , T* 

65  *( I )a2,*yMA*-y ( 1 1 

TAIL  LiMF(y,y,T*,i,o,oi 
no  66  1*1  , T* 

66  v(T)*?.*VMAy«Yf  II 

r*LL  LT5F(y, V, TX, 1 ,0,0) 

67  nn  6*  1*1,1* 

65  y ( I )*?.**xAy»* f 1 1 

CALL  LINF(y,y, tx, 1 ,0.01 
60  rnuTTMiE 
afTijo*) 

Fwn 

C 

«l)00nuTIME  BSOLV  (O.OA.OB.A^.ai  ,A01 
C CALCI/LATFS  BPfiTS  OF  OiJAOOATIC  FOIIATIOM 

c 

TF  (5P.FG.0.1  cnin  5 

BS(>A1  . SORTfAl**?  • 0,*A?*AU1)  /(2.AA21 
POINT  1 ,0,OA,09,A?,AI ,A0 

I FOOMAT  (OH  9 s,6F16.5) 

TF  (B.GE.Bi.iKd.B.LF.om  GOln  ? 

9*(-A1  « SOOT(At**P  • a.*A?*AU))  /(2,*A9) 

POINT  1 ,0,0A,PB,A?,Al,An 
? CONTINIIF 

PETIJON  ,, 


Pass  is  rc-ct 


% O««*0/»1 

OBtNT  1 ,B,ei , BA. t 7, i 1 . 10 
9FTlJO^ 

Fiin 

^U«BriuTT>'ie  )(Pii'TYf'<,>,7,TK) 

B»I^T8  CnriPDTMTES  OF  FQim.  VALUF  LI’^ES 

niMEMSinw  y ( 1 5 , y ( 1 1 
ORT^IT  8 
PPTNT  U,  z 
0BT*>iT  ?,  ly 
PPTNT  S 

optmt  1,  cx f K) , v(y ) ,K«i , ryJ 

1 P0P^‘lTflW  ,?FiP.T,^(?y,?Pt5.Tn 

2 rOP»*lTf1H  .TA.TW  POrVTT? 

3 POPMlTflw  ) 

a F0P“lTfQw  VIU.'F  a .PP,3) 

p pnPMxTftM  , 1 NurnnpniMiTE  pubst 

PETUPN 

PNT 

SUPPnuTINF  iPft,«,T,yw,yv,y,v«TA,ji,ty,NT,TI,vT 
PFAPrMFS  PnTTpM  ST^'E  he  RECTa'-RLE  FGP  fitVFN  VALtlE 

f>I«£K'STnM  yv  f 1 ) , Yv  f 1 ) , y f 1 T , y f 1 ) , vft  oo,?oo) 

TF  CA.r.T.T.ANp.H.tT.T.DB.A.LT.T.ANO.B.r.T.TT  1,7 

1 oi«yv(lAT 
i7t»XV{lA*l  ) 

*y»(Pl»(3i  )/?, 

T1»A 

TJlP 

TF  (TA.FO.n  7,3 

2 Pl«yv(TA*21 
T5«Vf IA*2, JA) 

CALL  8nL33f yy , T,P1 ,nt,Pl,Tl.Ti,T3) 
sn  Tn  t, 

3 TF  (TA.Eo.rn  «,s 

<1  PlBXV(TA«tT 

T3pvf r*»i . ja) 

CALL  SnL33fyy,T,Bl ,01 ,B1 ,T1 .T<,T3) 

on  TO  A 

5 TF  (V{TA-1  , JAl.LF.O.)  (JnTn  2 
TF  (V(TA*?,  JA)  .lf.o.t  P.nrn  « 

Bl»y V( T A-1 1 

Sl»yvfTA*2T 

T3iV(TA.1 ,JAT 
TU«vriA»?,jAT 

CALL  SnLUUfyy.T.Pl ,Ql,Bl,81»Tl,T7,T3,Ta) 
b MT»1 
Ty*iy*i 
y ( rxT«yy 

v(Tx)*yvrjA) 

TAsJA.I 
7 peTupn 
FNP 

SUPPOUTINE  MCfP,C,T,yV,yV,y»y*TA,JA,TX,NT,JJ,Vl 
8E*PCHFS  P.H,  STHE  ijr  PfCTANITLF  FHP  given  value 
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1 

f 


\ 


f 


> I 

I 

[ 


oi-FkSTOk'  yvfO.vvfn.vtn.’Ml.vfiPO.Pnn) 

TF  t«.RT,T.4vr.r.LT.T.np.u.(.T.T.4Kr».r.r,T.T)  1,7 

» 01»Yvf t4) 

01  »v  V ( t »♦  1 1 

YX« ( PI *01 ) /?. 

T1  «B 

T?»C 

TF  (J4.F0.n  ?,5 
? P1»VV(J*7?1 

T J»V  fI4*1  , Jk*y) 

C4LL  SniT5fXX,T,Ol,T1,P1,T1,T<:,TT) 

(jn  TP  B 

3 TF  «.*! 

u Pt»vv( J»»1 ) 

TSaVf I*4l ,J4-1 ) 

C*LL  S"L33fXT,T,Pl,Q1,Pl,T1.T^,T'«) 

(jri  TP  A 

F TF  (V(IA*1  ,J4-(  ).l.F.ft.1  r-OT'i  (? 

TF  f V(  T4*1  , J4  + ?1  .(.P  ,p,  T r,nTi-  u 
PlsYVt J«-1 ) 

F1«YV(.I4*?5 
T3»v ( I »♦! , J4-1  ) 

TU»Vf T**1 

PALI.  S^Laa  f X'',  T,P|  ,0t,oi,«l»T),T?,T3,Ta) 
fr  'JTa? 

T*»I**1 

Yf TXTxYVf T»*'  ) 

V ( T X 1 «** 

T4»I**1 
7 peT'.ipxj 
PVP 

c 

subspiittm?  cPfftf',T,yv,Yv,y.ViT4t.u.Tx,»T,ii,vT 
C FFAPPHFR  Tno  rTPF  np  PFfT^MULt  Fnp  tSTvFX'  W*L"F 

C 

PI-fMSTO*^  yvfO,vvfn,xtnf’fn,v(l0n,?00) 

TF  (r,r,T,T.»Nr.n.LT.T.np,p,LT.T,iKn,n,RT.n  1,7 
1 Pl«XV(T»1 
01 *xv  f T *♦ 1 1 
yx*(Pi*Oi )/?, 

T1*0 

T?«C 

TF  (lA.FP.tl  ?,3 
? P1»XV(T*+?1 

T3«Vf I4*?,J»71 1 

r*LL  snL33fxx,T,Pi  ,oi,pi,ti.t<!,t3i 
r,n  TP  Fi 

3 TF  (T4.F0.TT1  <I,F 
<1  PlaXV(Ti>n 

T3*vrr4-»,j4*t ) 

C*l-L  SFiL33f  XY.T.Oi  ,01,Pl,T1.Ti,T3) 
on  TP  F> 

5 TF  (V(  I*-l  , J»*n  .LF.P.)  pPT"  t 
TF  (V(T**?,J4*n.LF.fl.1  OPTI'  u 
PlaXVfT*-] 1 
FlaXv(T*7?T 
T3»V(I*«1 , Ji*l  1 
TOBY  f 1 1 

r*LL  soLauf yy.T.oi ,qi,oi,si»ti,t?,t3,tui 


I 


« ( T K 1 «y  It 
V(  Tx1ivv( J4*1  1 

PETU'PNi 

P»jn 


PuasnuTT^  E ''*t^,*,T,*v,yv,x.y#I4,J*,Ix,^T,,IJ,v5 

SP4PC«FS  L.h.  Flop  np  pp  c T » '"•P.l.  F POP  5IVFiy  Vil.IlF 

'iI«6MRTnt>i  yvfn,vvMi,y(n»'^(n,v(ioo,3nn) 

TF  (''.RT.T.*Mrs.».LT.T.np.n.LT.T.iX'n.4.r,T.T)  1,7 

1 oi*YV(J4) 

0 1 ■ y V t J 4 ♦ n 

yya{P1 *01) /?. 

T1  «4 
T?»r) 

TF  (j4.Fn.,)  p,i 

2 Pi«yv(.J4  + ?1 
T5lV(14,j4y?1 

r4LL  SnLXityy.T.Pl.Ol.Bl.TItTP.TH) 
r,n  TO  ^ 

1 TF  (J4.Fn.,IJ1  a,F 

u Plryvc  14-1 5 

73»yf  14, J4-1  1 

C4LL  SnL3ifyy,T,Bl,'51,Pl,T1»T<;,T3) 
tin  TO  t, 

5 TF  cv(T4,  j4-n  .lf.ti.t  r.PTp  g 
TF  fVtl4,J4»?l .LF.n.l  anio  u 
Bl»yv(j4-1 ) 

Rl«y V(  14*?1 
T5«vf  I4,>,T4-1  1 
T«»Vf  I4,.T4*?1 

r4LL  8'’L««fyY.T,t>l,'51,B1,S1/Tl,T?,T3,TU) 

6 MT«a 
Ty«!y«i 
y(TX3«yy(i4) 
yf Txiixx 
T49I4.1 

T PFTUO'g 

PNO 

8llP9nuTTVF  sn|33fxx,T,Bl,01,B),T|,T?,T3) 

pnLVtgnMTFL  ixirFBonixTiON  kITm  ^ ►■onFS 

xx9xy.Pi 

yt»o, 

y?90l -PI 

y39Ri-Pi 

009X2 

XX29X2*X2 

xxX9X3*y3 

TTpT.TI 

TT29T2-T1 

TT39T3-T1 

0T9xx2*x3-yx3*y2 

449(TT2*X3»TTx«X2)/nT 

q99(XX2*TT3-Xx3«7T2)/nT 

TX9.TT 

TZ9TT2-TT  65 


o c% 


•BIS  PlOE  IS  BIST  QUALin  PIUCT1CA0I»B 

Juan  ddc 


1 TV»(i4*y»>BH1*»K-TT 

TU»T**TY 
TP  (TU5 
? y^aXY 
TZaTV 
cn  TO  (1 
3 XtaXX 
T*«Tv 

<1  yii«(yiYX?W?. 

ll»  ( f /PQI 

XXsXd 

TP  (II. I T..nnB}  «i,i 
S xxaxx*Pi 
PpTllPKi 

PIwO 

^UB»nijTI*gF  ST'lO«f*X,T,Pl,Q1,Pl,Sl,Tt,T?,T3,Ta) 
pnt.yvnMfil.  T»uTFPonL*Tio^  *iIm  « Mnr>ES 

XXaXX.P) 

XlaO. 
yJaOt .P 1 
y3aP< -Pt 
XUaS I t 
PQax? 
yx?av?*y7 

yx3ax3*x3 

yx<iaya*xd 
yxv2ayx7*x7 
yxy3axx3*x3 
yxyuaxyu*xu 
TTaT-Tl 
TT?aT?.T1 
TT3aT3.T1 
TTUaTU.TI 

nTaOT3Tyxx7,vv7,y7,xyx3,yxX>xi,xxxu>yy<i,v<i) 
«4aOT3fTT?,»x?,x7,TT3,*y3,xi,lTu,yxu,xaWDT 
ppaOT3fxyx?,TT?,y?,«xx3,TT3»xl,xyxa,TTa,x«)/nT 
rcaOT3fiyx?,yy?,TT?,yxy3,xy3,iT3,xyxo,xx'i,TT«j/oT 
Txa-TT 
T?aTTJ-TT 

1 Tva ( ( * 4*y X tPP^ *x y*Pr T • xx-TT 
TllaTy*TV 
TP  (TiiT  P.B.x 
? y?axy 
TZaTY 
r,n  Tn  u 
3 xtaxx 
TXaTY 

u yi.ip(yi*x?l/?, 

'la4B3(  ( xx-XMi  /PQ) 
yxaXII 

TP  (".LT..nn‘;i  p,i 

3 yxaxv*P| 

9ETUPIY 
pmn 

c 

PU''CTTT'M  nTX(4»,9»,C),*?.«?»CP,43,e5,C3) 

C r4|.ril|.»  TPS  rPTFQ“TU4MT  DP  X BY  3 HXTBIX 


a.  3 I-J  pj, 

s>i4  y^f'i  r.«»..siir. 


bOtCHu 

Ol^ENSIUN  ►'*(lC0),H»(?0()1,«(u),rfg) 

INTECfcM  bUbttPPf  10),tJl’i(9P,10) 

OAT*  gL'^A/uM  / 

101  FOK-*T(tuF6.0) 

10?  A'JM***T(?OIa) 

110  AOS«**T(  ixl  ) 

1?1  fOBM*T(ln  , A 1 , iHHt ( , 1 }, 5h)  «, Ae . J, U ( Au , 6A, iHh* ( , 1 J, JH)  BfrA.i)) 

131  AurtMATll"  . A1  , 5>*«»  ( » 1 i.  JH)  »,A»‘,3,u(4a,o*,3>’r«r(,13.3«)  •.AT. 3)) 
?0?  A.jh«at(ih  , Ju.SfoA.SIu) ) 

HEAD  10^,^*,^r,^EP>«,^tP^ 

PEAO  10?,  10,  I^,  Jl’.,  JN 

NIM 1 BNXa 1 
S»(«1«SV»1 

«EA0  101  , (>tl  ( n , lal  ,^A  ) 

MEAD  10 1 . (»T ( J) , J«  I ,ST  ) 

pEaO  io?,(teo«(i,j),j«i,NtPx).i«i,NrMi) 

«EAo  io?,((bo?(i,j),j«i,Ntp?i,i«i,N«"n 
HXO*HA ( 10} 

UO  100  lain, IN 

HX ( l»i-lu)aHX ( I ) 

00  lao  I«ai,NEP? 

IP  (a02(I,lK).Ll.?*Ju)  bOZ( t, lK)a?*JO 
If  »0?(I,  J-<)a?*J^ 

10  0 UO?  ( lAl>IO.  IK)abl)Z(  I , 1«  )•?•  ( JU«1  ) 

NXalN*l»IO 
NXMlaNXal 
MYOaHYl JO) 

00  laO  JaJu.JN 
HV(jAl«JU)aHY(J)-HYI' 

00  160  JKai.^EPK 

lA  (aDM(J,J«).LT.?»Iu)  bO«(  J,  JMa?*IO 
lA  (OOXC  J,  JT)  .i,T  ,?*1M  bO'YCJ,  JK)a?*lN 
160  aOH(jAt>J(),  JK)aT>)P(  J,  JK)«?*(  iu>n 
ixYa  JN4 1 • Jo 
NYMlBMY-l 
ppiNi  no 
«laO 

•xCouaNXTl/SO 

6oro  ??? 

??0  NCOLa^COL-l 
???  IN»asCOL*50 

??S  TlaMi*! 

T?aai«lNA 

PPI**}  1?1 , (hcmk.h.hi  (X  ) .KaKl  ,K?,S0) 

lA  (Tl.Cb.90.Ort.Al.EW.NX)  bOTu  ??P 
iA  (A?. CO. NX)  ??0,??9 
??9  PPINT  UO 
AlaO 

NCOLaNYrtl/90 
SOTO  ?3? 

?30  NCOLaNCOft 


I 


I 

t 

I 


i 

I 


( 


1 

I 


2J2  lMi»NCUt*SO 
2JS  ■(li'Htl 
i^2iki*Imr 

PR  INI  1JI,{RLM«,K,"TI»),**M|»2,50) 

IF  (KJ  .kUI.50.LiR.Kl  .EW.NV  ) tiUTO  2J«> 

IF  (k2,EQ,nv)  2Jo,2JS 

2j<j  PRIM  no 

OC’  eso 

250  PRINT  202.J#(ttOR(J.K),»»l,s£PR) 

PRINT  MO 
OC  260  J«l,N«Ml 

260  PRINT  202.  J,  1602(  J.x  ) ,R«MNtP2) 

C*LL  PLOTS  (M.0,90) 

C*LL  PLOT  (0,,-ia.,-l) 

C*Lt  P|.l)T(0.,0.5,-4) 

XM/|XSHIt(M).MX(  1 ) 

YHAXBHVINYJ.hYM) 

F IRSTX6HX( 1 ) 

OCLT*x«*P*«/lO, 

FIB5TY*HY(n 

OeLT*Y«UEl.f*« 

Y*xaYRix*10,/XM6x 

C*LL  *xIS(0., lO,,6H2-*xIS,*b,Y*«.0.,FiRSTY,UEUT*Y) 

CAUL  AXISTO., t0.,6HW.»xlS,-6, I0.,27  0.,FIRST  x.OELTAX) 

xCJIaFIRSTx 

X(«)aDELTAX 

Y(3)aFIRSTY 

Y(P)ai3ELTAY 

00  500  lal.NYPl 

YMIaHYdYl) 

y(2)aHY(n 
UO  UUO  Jal.NbPR 

IF  {t»0Rll,J).CT,2xNx)  Goto  500 

LattOR(I,J)/2 

KaSORC  I. Jl/2.  Y .5  • L 

X ( 1 )aXMXx>HX (L) 

IF  (a.gE.M  goto  350 
Fr2)ax(  1 ) 

GOTO  390 

350  X(2)aXP*x-nx(LYl ) 

390  call  LlNE(Y.x,2.1.0.b) 

900  IF  (S0R(I.J).GE.2*NX)  C(>TU  500 
500  CONTINUE 

00  700  lalfNXMl 
X ( naxPXX>MX  ( I ) 

X(2)aXHAX>HX(lYl) 

DO  600  Jal.hEPZ 

IF  (60Z(I*J).GT.2*NY)  goto  700 
LaB0Z(i«J)/2 
Ka60Z(l.J)/2.  Y .5  « L 
Y(t )ahY(L) 

IF  (K.GC.l)  GOTO  550 
V(2)aY(l) 

GOTO  590 

550  TtZIaMYtLYl) 
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I 
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590  call  LINe(YiX,2.1,0,0) 

600  IF  (dOZ(I«J).GE.2*‘YY)  GOTO  700 
700  CUNTiNUt 

call  PLOT(tZ.0,0.0,«0) 

STOP 

iNO 


PAOGRAM  EPSCHK 

dimension  B{I»9,  J'J'J)  ,Z(9<»,  199)  ,EP8a(  >69,  tO),tPSZ(99, 10) 
iNTtGEP  80B(199,10),8DZ(99,IO) 

DATA  EP SR/ 199U*0.U/, I PS 7/990*0.0/ 

DATA  «/t970l*-l.O/,Z/1970l*»1.0/ 

102  FCRmaT(201u) 
lOa  FOR"ATUeF5,0) 

150  format  (lni,)x,52MN  n*  NT  i«(n*,nv)  2(*<x,ny)) 

IS5  FORMAT  (IH  ,SlS,2Fie,j) 

161  FORMAT  (Iho,16«,2«nu) 

162  format  (|M0,12»,10Hi)ieLtCT«lC) 

16J  format  (ImO,  to*,  ISmTnC'JNSISTInCIES) 

175  FORMATtlMl,5x,lMA,2<i(I«,l«)) 

176  F0HMAT(6*,1m*,2«(1«,1«)) 

160  F0RmaT(6*, IH* ,2«(Sm..-a.)  ) 

165  FOHmaT(1*,Iu,2m  ♦,2<*F5,1) 

READ  t02,N*,Nr,NEPR,NfcRZ 

NXM 1 6N*a 1 
NVMlaSY-1 

mEAO  102, ( (60R( I • J) , Jsl .NEPR},I61,6YM|) 

2 READ  102,><U,hn,ka 

READ  I on,  (tPSR(«li,*  ) ,«*  1 , *K1 
IF  (AO.EG.XN)  tUTy  « 

R I 6*0* 1 
DO  J 1**1, RN 
DU  3 R*1,RR 

i E*SR(l,R)*eRS»(RU,R) 
a IF  (rr.EO.nEPR)  GOTm  5 

RRaRR*  1 
00  *5  l*RO,RN 
DO  «5  J*RP,nCPm 
U5  EPSR( I, J)*EP6R(  1,RR) 

5 IF  (RN.LT.nymI)  COT''  2 

READ  102,((RDZ(1, J),J*1,«EP7),1*1,nxm|) 

6 RE  AO  1 02,  RLI,  RN,R* 

READ  lOn, (EP8Z(RU,R),«*1 ,«Rl 
IF  (RO.EO.RN)  GOTO  6 
KlaROtl 
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II  .1  » g* 
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‘ ro  OOC 
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00  7 I«Hl,M\ 

00  7 Kai««K 

7 tPSZd.KjtfePSZfKU,-^) 
tt  IF  iRK.EW.NfcPZ)  GOTti  << 

KPsKKt  t 
00  89  I»«0,«N 
t)J  89  J»nP,N.EPZ 
89  EPSZ ( 1 • J)aEPSZ ( 1 >kk  ) 

9 IF  (kn.lT  .NXHi  ) gOTii  t 

00  ao  lal.sMrtj 

HaO , 

KjaO 

00  ao  Jal,NfcP>* 

IF  (M-^a.GT.u.)  *aa«a*» 

KJana*! 

M a0OP(I,JJ/a.-l. 

Hfan 

IF  (Ka.LT.Ki)  goth  a.)  j 

00  10  Kaal,«a 
10  ><(R.I)afcPSPlI.J) 

ao  continue  5 

00  «0  Ial,N*8i  ; ; 

"ao,  ^ 

8aaO  " 

DO  «0  Jal.NEPZ  j 

IF  ih-nz.gt.o, j 

KlaKitl 

" ■8DZ(I» J)/a,»J  . 
aaa" 

IF  ti'a.CT.Kl)  UOT'J  ao 
00  50  KiKifKa 
50  Z(I«8)aePSZll#J) 

UO  CONTINUE 
NaO 

PHINT  150 
00  60  lal,NX") 

00  60  Jal.NTMl 

IF  IP(I#J).Fw.Z(I,J)T  -JTO  60 

NaNf  1 

ppiNT  iss,N,i,j,-(i,j),afi, j) 

IF  (N.St.lOOO)  SJT(j  7 0 

60  continue 

IF  (N.DT.O)  UuTU  70 
PPINT  161 
PPlNt  16a 
PPINT  165 
70  continue 
210  "2aO 
"1650 

220  "26"2«50 

IF  ("2.tE,Nt)  SOTO  225 
"16NY»50«"2 
"26NV 
225  N260 
250  NlaN2tl 
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N2aN2t2U 

IF  (N2,ST.mx)  n2».m» 

PrtINT  175» 

PRINT  leo 
OJ  2b0  Jal.**! 

I»P2*l-J 

250  PRINT  195. 1.  («(■'.  I ) jXaNl  ,N2) 
PRINT  lao 

PRINT  176,  (K,i»»ni,n2) 

IP  CN2.CT,nx)  goto  230 
IF  («2.LT.ny)  goto  220 
260  M2a0 
Ml  a50 

270  M2aM2t5n 

IF  (m2,lE,n>j  goto  275 

mi»nv»5o-m2 

MgmH* 

275  N2a0 
200  NlaN2tl 

IF  (n2,GT,nx)  N^a^x 
PRINT  175,(x.NaAii,N2) 

PRINT  190 
DU  300  Jal.Mi 
IaM2tl»J 

300  PRINT  165. 1. (2(N. n .NaNl .n2) 
PRINT  160 

PRINT  176. (x,KaNi ,N21 
IF  IN2,LT,nx)  goto  260 
IF  (M2.LT.Ny)  GOTO  270 
STOP 

EM) 
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